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ABSTRACT 


Ccmtributions  of  the  anode  to  Magnetopiasmadynamic  (MPD)  thruster  performance 
are  ocmsidered.  energy  losses  at  this  electrode,  surface  erosion,  and  sheath/ionization 
effects  must  be  controlled  in  designs  of  practical  interest.  Current  constriction  or  ^tting 
at  the  anode,  evolving  into  localized  surfece  damage  and  considerable  throat  erosion,  is 
shown  to  be  related  to  the  electron  temperature’s  (T.)  rise  above  the  gas  temperature  (T,). 
An  elementary  onenlimensional  ctescription  of  a  coUisional  sheath  which  highlights  the  role 
of  T,  is  presented.  Cmnputations  to  model  the  one>dimensional  sheath  are  attempted  using 
a  set  of  five  coupled  first-order,  nonlinear  differimtial  equations  describing  the  electric  field, 
as  well  as  the  species  current  and  number  densities.  For  a  large  temperature  nonequilibrium 
(i.e.,T,>  >To),  the  one-dimensional  approach  fails  to  give  reasonable  answers  and  a 
multidimensional  description  is  deemed  necesssary.  Thus,  anode  spotting  may  be 
precipitated  by  the  elevation  of  T,  among  other  factors.  A  review  of  transpiration  cooling 
as  a  means  of  recouping  some  anode  power  is  included.  Active  anode  cooling  via 
transpiration  cooling  would  result  in  ( 1 )  quenching  T„  (2)  adding  "hot"  propellant  to  exhaust, 
and  (3)  reducing  the  local  electron  Hall  parameter.  However,  significant  technical  problems 
remain. 
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I.  INTRODUCTION 


Several  types  of  space  flight  propulsion  systems  have  been  developed  over  the 
years.  These  include  chemical,  nuclear,  electric  and  solar  propulsion.  The  majority 
of  space  thrusters  to  date  have  been  chemiod  rockets.  Although  the  Qiinese  used 
rockets  over  800  years  ago,  true  development  of  rocket  fH'opulsion  took  place  during 
this  century  [Ref.  1].  Chemical  thrusters  give  high  thrust-to-weight  ratios,  larger  than 
unity,  and  have  been  fully  developed  in  the  form  of  space  launch  vehicles  and 
attitude  control  thrusters.  In  contrast,  other  propulsion  systems  have  been  developed 
only  to  the  proof-of-ooncept  stage,  and  essentially  remain  at  this  stage  of 
development.  Nudear  propulsion  was  studied  with  the  NERVA  (Nuclear  Engine  for 
Rocket  Vehicle  ^plication)  thruster  in  the  1960’s,  and  abandoned  [Ref.  2:pp.  517- 
519].  Electric  propulsion  flights  during  the  1960’s  included  the  U.S.  SERT-1  (Space 
Electric  Rocket  Test)  in  1964  and  the  U.S.S.R.  Yantari-1  rocket  in  1966.  Solar- 
electric  propulsion  was  demonstrated  via  the  SERT-2  rocket  in  1970,  powering  the 
electric  thruster  from  power  generated  by  solar  cells.  Further  electric  propulsion 
research  flights  in  the  1980’s  induded  the  U.S.  Navy’s  NOVA-1  satellite  in  1981,  and 
Jsqpan’s  MS-T4  satellite, launched  from  the  Space  Shuttle.  Beyond  this,  nonchemical 
thrusters  have  only  been  used  in  auxiliary  roles,  such  as  station-keeping  and  attitude 
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control  on  geo^chronous  satellites.  NASA’s  Project  PATHFINDER  in  the  mid- 
1980*s  proposed  the  use  of  a  megawatt-level  electric  plasma  thruster  for  a  manned 
Mars  mission.  However,  development  of  this  project  was  never  funded. 

In  comparing  the  different  propulsion  schemes,  a  primary  performance  indicator 
is  specific  impulse,  defined  as  the  ratio  of  thrust  to  the  rate  of  propellant  usage,  or 
alternately,  propellant  effective  exhaust  velocity  (u.),  divided  by  the  sea-level 
gravitational  constant,  (g^). 


I 


•p 


So 


sec 


(1) 


Chemical  rockets  are  inherently  limited  in  performance  by  the  total  energy  available 
in  the  fuel/oxidizer  combustion  process,  so  that  the  total  enthalpy  available  for 
conversion  into  exhaust  kinetic  energy  is  limited.  Exhaust  velocity  is  also  limited  by 
material  heating  limitations  of  the  combustion  chamber  and  nozzle  throat,  and 
"frozen  flow  Losses"  (unrecoverable  energy  deposition  in  internal  modes  of  the  gas) 
[Ref.  3;pp.  4-5].  Peak  specific  impulse  for  liquid  chemical  propellants  is  presently  on 
the  order  of  450  seconds.  This  capability  is  completely  sufficient  for  the  tasks  of 
launch  to  low  earth  orbit  (LEO),  attitude  control,  station  keeping,  and  such  missions. 
However,  for  missions  such  as  manned  interplanetary  exploration,  chemical 
propulsion  can  be  shown  to  be  clearly  inadequate.  A  comparative  analysis  of  a  Mars 
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Biuakm  using  chemical  and  electric  propubion  systems  shows  the  large  difference  in 
mass  payload  ratio  (final  mass/initial  launch  mass)  for  the  two  systems.  A  chemical 
system  using  a  hi^  imimise  Hohman  ellipse  trajectory  delivers  a  maximum  of 
afi^ifoximatefy  10%  to  18%  of  the  launch  mass  to  the  Martian  surface  [Ref.  4:p.  US]. 
In  conqKurison,  an  electric  system  using  a  low  impulse  spiral  trajectory  could  deliver 
from  20%  to  60%  of  the  launch  mass,  depending  on  the  desired  transit  time.  Each 
mission  assumes  transit  from  low  Earth  orbit  to  Mars  orbit.  An  electric  propulsion 
^tem  would  still  need  a  high  thrust  propulsion  system  to  reach  the  Martian  surface 
Ref.  S:pp.  344-346].  The  large  difference  in  payload  ratio  is  due  to  the  much  larger 
exhaust  velocity  and  more  efficient  use  of  fuel  by  electric  propulsion.  Thus,  some 
form  of  electric  or  hybrid  electric  thruster  would  seem  to  be  in  order  for  such 
interplanetary  missions.  However,  due  to  the  low  thrust-to-weight  ratio  of  electric 
thrusters,  they  must  be  launched  into  orbit  by  other  means.  Their  usefulness  is 
restricted  to  space  thrusters,  not  to  launch  systems. 

With  specific  impulses  of  as  high  as  10,000  seconds,  elearic  propulsion  offers 
the  performance  envelope  needed  for  manned  interplanetary  missions.  Electric 
propulsion  is  divided  into  three  types  of  thrusters:  electrostatic,  electrothermal,  and 
electromagnetic.  The  type  relevant  to  this  work  is  the  magnetoplasmadynamic  (MPD) 
thruster,  an  electromagnetic  propulsion  tystem  that  utilizes  the  Lx)rentz  force  created 
by  an  electric  current  together  with  its  induced  magnetic  field  to  propel  a  gas  that 
has  been  heated  to  the  plasma  state.  According  to  electromagnetic  theory,  a 
conductor  carrying  a  current  produces  an  induced  magnetic  force  perpendicular  to 
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the  current.  The  applied  electric  field  and  its  induced  magnetic  field  interact  to 
produce  the  Lorentz  force  if  •  JxS)  perpendicular  to  both  fields  on  the 
conductors.  This  briefly  summarizes  the  concept  behind  the  "self-field"  MPD 
accelerator  [Ref.  2:pp.  485-486].  MPD  performance  is  enhanced  by  adding  magnetic 
coils  to  the  thruster,  thus  strengthening  the  magnetic  field  and,  as  a  a>nsequence,  the 
Lorentz  force  and  thrust.  This  thruster  is  appropriately  called  an  "applied-field"  MPD 
thruster.  MPD  thrusters  have  shown  specific  impulses  of  up  to  7,000  seconds  and 
efficiencies  as  high  as  70%  [Ref.  6:pp.  2-3].  Performance  of  MPD  thrusters  is  limited 
by  several  factors,  including  electrode  erosion,  current  spotting,  frozen  flow  losses, 
and  electrode  power  deposition.  Specifically,  anode  power  deposition  is  the  single 
largest  power  loss  mechanism  in  MPD  thrusters  operating  at  submegawatt  power 
levels  [Ref.  7].  In  the  following  work,  we  review  and  analytically  model  the  MPD 
anode,  including  the  sheath  and  anode  potential  drop. 
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II.  LITERATURE  REVIEW 


Anode  losses  significantly  limit  magnetoplasmadynamic  (MPD)  thruster 
performance.  Much  effort  has  been  placed  on  characterizing  these  losses  and  on  the 
nature  of  power  deposition  in  the  anode  (Refs.  8-14].  As  much  as  80%  of  thruster 
total  power  may  end  up  being  deposited  in  the  anode  at  sub-megawatt  power  levels 
[Refs.  8,15].  This  power  deposition  together  with  current  constriction  at  the  anode 
surface  present  serious  problems  to  thruster  cooling  and  performance,  as  well  as  to 
anode  lifetime.  Before  any  practical  design  can  be  achieved,  a  more  thorough 
understanding  of  the  phenomena  at  the  anode,  particularly  the  anode  sheath,  must 
be  gained.  Studies  have  shown  that  the  anode  power  fraction  depends  on  thruster 
power,  current,  mass  flow  rate,  and  the  parameter  J^/ih  [Refs.  8,12,13,16).  It  has  also 
been  shown  that  the  anode  fall  voltage  is  inversely  proportional  to  anode  current 
density  [Refs.  13,16].  It  is  believed  that  a  better  understanding  of  the  role  of  an 
elevated  electron  temperature,  of  current  flow  dimensionality,  and  of  current 
unsteadiness  are  prerequisites  for  the  evolution  of  any  practical  MPD  thruster  design. 

Computer  codes  that  accurately  describe  observed  data  from  steady-state  MPD 
thrusters  have  been  developed  [Refs.  17-19].  However,  these  codes  do  not  adequately 
describe  observed  data  from  quasi-steady  thruster  experiments.  It  has  been  suggested 
that  the  lack  of  proper  electrode  modelling  (i.e.,  sheaths  and  fall  potentials)  in  these 
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codes  may  explain  this  discrepancy  [Ref.  6].  Limited  analytical  work  has  been  done 
in  modelling  the  sheath  and  ambipolar  regions  at  the  anode,  influenced  perhaps  by 
the  difficult  set  of  coupled,  nonlinear  partial  differential  equations  involved.  Hugel 
[Ref.  12]  and  Subramaniam  [Ref.  20]  address  the  influence  of  the  sheath  region,  but 
do  not  model  the  electric  field,  temperature,  or  sheath  fall  voltage. 

Given  the  minuscule  extent  of  the  sheath  versus  thruster  anode  curvature,  the 
problem  at  first  appears  one-dimensional  in  nature.  A  one-dimensional,  collisional, 
equilibrium  solution  can  satisfactorily  reproduce  the  observed  electric  field  and 
charge  density  distributions  for  the  entire  sheath  and  ambipolar  regions  for  a  sheath 
where  the  electron  temperature  equals  that  of  the  heavy  species  [Ref.  21].  However, 
this  model  cannot  describe  any  decrease  in  current  density  away  from  the  surface,  or 
current  constriction,  at  the  anode  surface  which  might  be  necessary  in 
nonequilibrium.  A  two-dimensional  model,  developed  by  Biblarz  and  Dolson  [Ref. 
14],  represents  these  phenomena  and  predicts  the  voltage  drop  in  the  region.  It  is 
shown  that  the  sheath  must  account  for  a  majority  of  the  anode  voltage  drop,  and 
that  the  sheath  extent  must  be  greater  than  the  Debye  length  [Refs.  14,21).  Thus, 
a  combination  of  one-  and  two-dimensional  approaches  appears  to  better  describe 
sheath  behavior.  Incorporation  of  modelling  of  this  sort  may  improve  the  ability  of 
the  computer  codes  cited  above  to  properly  describe  quasi-steady  thrusters. 

Next,  a  description  of  the  anode  region  is  presented  in  order  to  delineate  some 
of  the  possible  effects  of  temperature. 
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III.  Am>DE  DESCRIPTION 


A.  THRUSTER  GEOMETRY  DESCRIPTION 

The  majority  of  plasma  thrusters  to  date  have  consisted  of  a  central  cathode 
rod  surrounded  by  an  annular  shell  anode,  as  shown  in  Figure  1  [Ref.  23].  The 
thruster  illustrated  is  sufficient  to  produce  needed  thrust  at  current  levels  above  one 
kiloamp.  Below  this  level,  an  external  magnetic  field  produced  by  an  annular  magnet 
is  needed  to  ensure  sufficient  Lorentz  force  on  the  plasma  propellant  to  meet  thrust 
requirements.  [Ref.  8].  As  illustrated  in  Figure  1,  the  jxB  body  force  simplifies  into 
an  axial  (jgB^)  body  force,  whidi  provides  direct  electromagnetic  thrust  ("blowing"), 
and  a  radial  body  force,  which  provides  electromagnetic  compression  of  the 

plasma  and  a  subsequent  pressure  force  along  the  cathode  surface  ("pumping”).  [Ref. 

6] 

A  notable  exception  to  this  geometry  is  the  Stationary  Plasma  Thruster  (SPT), 
a  design  firom  the  former  Soviet  Union.  The  SPT  is  an  example  of  a  plasma 
propulsion  system  known  as  a  Hall  Current  Plasma  accelerator.  An  electric  field  is 
applied  axially  to  a  stream  of  flowing  plasma,  in  addition  to  a  magnetic  field  with  a 
strong  radial  component,  which  is  applied  by  an  external  electromagnet.  When  the 
axial  electric  field  is  applied  and  a  current  flows  through  the  plasma,  an  azimuthal 
component  of  current  is  induced,  i.e.,  the  "Hall"  current. 


7 


PROPtLLANT 


Figure  1  >  Magnetopiastnadynamic  ^ . 

on  Plasma  Indicated.  (Ref.  23] 


Thrust  is  imxtuced  electrostatically  accelerating  the  ions  in  the  plasma,  as  well  as 
through  the  induced  Lorentz  force  mentioned  above.  A  strong  radial  magnetic  field 
is  lulled  to  the  plasma,  whose  properties  are  controlled  to  make  the  electron 
Lammr  radius’  small  compared  to  the  mean  free  path’,  while  the  ion  Larmor  radius 
is  comparatively  large.  As  a  consequence,  electron  mobility  in  the  axial  direction  is 
greatly  reduced.  Thus,  the  electric  field  energy  is  given  mainly  to  the  ions,  pr  ng 
axial  ion  acceleration.  Collisions  with  neutral  particles  serve  to  accelerate  the  entire 
neutral  plasma.  [Ref.  24] 

A  pair  of  the  final  prototype  design  developed,  the  SPT-100,  have  been 
acquired  NASA  recently  from  Fakel  Enterprise  in  Kaliningrad,  Russia,  and  are 
undergoing  performance  evaluation  at  the  Jet  Propulsion  Laboratory.  Designed  at 
the  Kurchatov  Institute  of  Atomic  Energy  (LAE)  in  Moscow,  USSR  in  the  1960’s, 
smaller  versions  of  the  SPT-100  (SPT-50  and  SPT-70)  were  flown  beginning  in  1972*. 
A  specific  impulse  of  1,600  seconds  and  S0%  efficienqr,  as  well  as  space  flights  of 
fifty  similar  thrusters  is  claimed.  The  specific  operational  characteristics  of  the 
thruster  are  not  well  understood  presently.  Bohm  diffusion  of  electrons  and  a 
phenomenon  called  ”near-wall  conductivity"  have  been  proposed  to  explain  the 
thruster’s  operation.  This  thruster  is  shown  in  Figure  2.  [Ref.  25] 

‘  Larmor  radius  is  the  radius  of  the  helix  traversed  by  a  charged  particle  moving 
in  a  magnetic  field. 

’  Mean  free  path  is  the  distance  traveled  by  a  particle  before  making  a  collision. 

The  suffix  (i.e.,"-70")  indicates  the  characteristic  diameter  of  the  thruster  in 
millimeters. 
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B.  ELEMENTARY  SHEATH  FORMULAE  DESCRIPTION 


1.  Disaissimi 

Voltage  losses  and  anode  power  deposition  account  for  most  of  the 
inefficiency  of  plasma  thrusters.  In  order  to  understand  these  losses,  the  anode  region 
must  be  understood  and  related  phenomena  explained  and  modeled.  As  shown  in 
Figure  3,  a  substantial  drop  in  voltage  occurs  in  a  short  distance  from  the  anode 
surface. 


Cathode  Anode 

Figure  3  -  Electric  Field  Between  Two  Electrodes,  Including  "Positive  Column".  (Ref. 
27J 

The  anode  fall  region  may  be  divided  into  two  parts,  the  sheath  and  ambipolar 
regions.  The  plasma  attempts  to  adjust  itself  near  electrodes  so  as  to  shield  the  main 
body  of  the  plasma  from  the  electric  field  (Ref.  26].  The  sheath  is  the  region  closest 


II 


to  the  enode  surface  within  ^Ndikh  the  ion  and  electron  number  densities  are  unequal, 
with  the  electrons  dominating  the  region.  A  high  electron  sfMice  charge  exists  at  the 
anode  surface.  This  is  caused  by  the  anode  collecting  incoming  electrons  in 
convicting  the  arc  current  with  the  cathode.  Positive  ions  are  produced  within  the 
sheath  by  electron  impact  of  neutral  gas  molecules,  and  the  ions  are  repelled  toward 
the  cathode.  At  the  cathode  end  of  the  anode  drop  region,  the  density  of  positive 
ions  is  high  enough  to  almost  neutralize  the  electron  space  charge,  thus  forming  the 
positive  column  or  core  plasma.  The  essential  positive  ion  current  is  created  in  this 
way  near  the  anode.  A  more  complete  description  of  this  process  may  be  found  in 
Cobine  [Ref.  27]  and  von  Engel  [Ref.  28].  A  fundamental  characteristic  of  plasma 
behavior  is  its  tendency  toward  electrical  neutrality.  Whenever  local  charge 
concentrations  arise  or  external  potentials  are  introduced  into  a  system,  these  are 
shielded  out  in  a  distance  known  as  the  "Debye  length"*.  This  distance  must  be  much 
smaller  than  the  system  dimenuon  for  the  ionized  gas  to  be  considered  a  plasma 
[Ref.  29].  Equation  (2)  gives  the  Del^e  length  [Ref.  26]. 


niie  Debye  length  effectively  describes  the  radius  of  a  shell  around  a  charged 
particle  outside  of  which  the  potential  of  the  particle  is  not  seen. 
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Aaotiier  disumce  of  interest  is  the  electron  mean  free  path,  or  distance  traveled  by 
a  particle  before  making  a  collision.  Equation  (3)  is  from  a  derivation  of  Lin,  Resler, 
and  Kantrowitz  [Refs.  3031]  giving  the  mean  free  path,  with  being  the 
proximate  sheath  length. 


X 


0.12 


_ 1 _ 

n.(eV3kT)*ln(k.eV3kT) 


(3) 


Since  the  sheath  extends  at  most  a  few  mean-free  lengths  from  the  anode  surface, 
curvature  of  the  anode  does  not  affect  the  governing  equations  in  high  pressure 
discharges.  Thus,  the  region  may  be  described  in  one  dimension,  the  distance  V 
from  the  anode  surface.  While  the  Debye  length  is  sometimes  assumed  as  the  sheath 
extent.  Reference  22  showed  that  the  sheath  thickness  is  a  frinction  of  the  anode  fall 
voltage  and  the  electron  temperature.  Equation  (4)  gives  the  £q>propriate  form. 


e  • 

O  • 


\  eiL 


(4) 


An  example  case  with  a  fall  voltage  of  100  volts  gives  a  sheath  extent  of 
*  2 .352x10'^  m.  This  compares  to  a  computed  Debye  length  of 
Aj,  «  1. 690x10*^  m.  Therefore,  the  sheath  can  be  an  order  of  magnitude  larger 
than  the  Debye  length. 
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Nasser  [Ref.  32]  dismisses  an  elementary  theoretical  approach  to  the  glow 
d^diarge  problem.  He  suggests  a  set  of  four  one-dimensional  ordinary  differential 
equatiaos,  including  the  electron  and  ion  current  and  number  density  equations,  in 
addition  to  Poisson’s  equation.  Most  solution  attempts  have  failed,  with  the  boundary 
condititHK  being  identified  as  the  culprit  A  similar  attempt  for  the  plasma  thruster 
is  discussed  below. 

2.  Simplified  Formulation 

The  steady  probe  equations  are  first  written  [Ref.  21]  in  their  simplest 
form.  The  anode  is  assumed  to  operate  as  a  heavily  biased  probe,  which  is  true  for 
low  enough  currents  when  the  anode  is  not  a  source  of  ions.  Whenever  the 
temperature  can  be  considered  fixed,  the  energy  equations  are  implicitly  satisfied 
and,  since  ion  inertia  is  neglected,  the  resulting  set  consists  only  of  two  species 
continuity  equations  and  Poisson’s  equation.  These  equations  are  written  in  terms 
of  y,  which  is  the  coordinate  outward  from  the  planar  positive  surface.  Constants  and 
variables  are  listed  in  Table  1. 
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I  y.1 i  31  Oa<KTt  ajgii  ^  iiTTs 


a...diaracteristic  length  of  i^asma  n....species  number  density  at  core  plasma 

Di.,..4^cies  division  coefficient  N...total  number  density 

e...elementaiy  cimrge  constant  T...temperature 

E...electric  field  ...neutral  species  temperature 

Eq ...electric  field  at  anode  surface  a. ..two-body  recombination  coefficient 

E»...electric  field  at  core  plasma  e^,  ...permittivity  constant 

j  ^ •  ...species  current  density  v  ...ionization  coefficient 

J. ..total  current  |ij^,«...species  mobility  coefficient 

k...Boltzniann’s  constant  ^^^...anode  fall  potential 

K. ..current  parameter  X..jnean-firee  distance 

....spedes  number  density  A4...0ebye  length 

fi,...time  rate-of-change  of  n.  X.  ...Sheath  thickness 

Note:  Species  subscripts  denote  ions  (i)  and  electrons  (e). 
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dy 

(5) 

j.  -  eM.n.E  ♦ 

(6) 

dy  e. 

(7) 

J  -  j.  ♦  j. 

(8) 

(9) 

Here,  the  f  s  are  spedes  contributioiis  to  the  total  current  densi^.  The  existence  of 
negative  charges  as  free  electrons  is  [^otal  in  the  formulation.  Next,  the  Einstein 
relation,  etpiation  (9),  is  introduced  to  write  the  mdlMlities  in  terms  of  the  diffusion 
coefficients.  We  assume  that  the  diffiisicm  coefficients  remain  constant  in  the 
proUem. 

Equations  (5)  and  (6)  are  next  solved  fm  dnj.  Jdy.  The  species  current 
density  equations  are  found  from  the  net  reaction  rate  of  the  plasma.  Equations  ( 10) 
ami  (11)  combine  to  produce  space  derivatives  for  species  current  density. 


ft.  »  VB,  - 


an,n. 


(10) 


dy 


(11) 


Combiniiig  equations  (SHU)  produces  a  set  of  five  coupled,  non-linear  differential 
equations  describing  the  sheath.  These  are  nondimensionalized  to  adjust  all  variables 
to  the  first  order,  and  are  rewritten  below  as  equations  (12)-(16),  with 
nondimensionalized  variables  denoted  Nondimensionalization  can  be 

accomplished  as  follows:  The  species  number  densities  n«n„  are  divided  by  their 
values  at  infinity  to  produce  output  from  the  anode  surface  to  unity  at  the  ambipolar 
boundary.  The  current  densities  j«ji  are  divided  by  the  total  current,  allowing  the 
output  to  show  the  "mirror  behavior"  of  the  two  currents.  The  electric  field  is  divided 
by  the  initial  anode  value  to  give  output  starting  from  vmty  at  the  surface  and 
decreasii^  to  the  final  core  field  value.  The  variable  Y  is  divided  by  the 
characteristic  length’  of  the  plasma,  "a",  produdi^  These  corrections  allow  all 
output  to  vary  in  the  range  from  zero  to  one,  as  a  function  of  distance  from  the 
anode. 


niie  characteristic  length  is  defined  so  as  to  cancel  the  multiplying  factor  in  the 
electric  field  equation,  (14),  (a«  1.107  x  10^).  This  allows  a  physical  interpretation  of 
the  ion/electron  number  densities,  as  well  as  the  decay  rate  of  the  electric  field. 
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Attenqns  to  solve  tlus  equation  set  using  tlM  conqiuter  code  discussed  below  shows 
the  set  to  be  extreme^  sensitive  to  initial  a>nditiom.  The  computer  code  solver  uses 
a  "nuffching"  scheme  from  the  anode  to  the  undisturbed  plasma.  The  initial 
omditioiis  are  chosen  to  produce  the  electric  field  potential  drop  observed  in  actual 
thrusters.  First  and  second  space  derivatives  of  the  electric  field  are  used  as 
diagnostic  checks  to  ensure  reasonable  output  values  and  indicate  instability  of  the 
int^pation  process.  Figure  4  shows  the  required  resulting  curves  for  the  electric  field 
and  its  first  and  second  derivatives. 
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Ecker  characterizes  the  plasma  at  the  anode  as  a  double  sheath,  with  the  inner 
section  called  the  "inertia  sheath",  and  an  outer  section  called  the  "energy  loss 
section".  The  inner  section  shows  a  potential  rise  of  the  order  of  one  volt,  with  the 
outer  section  showing  the  ej^nential  potential  drop  shown  in  Figure  4.  While  this 
double  sheath  may  in  fact  describe  the  actual  sheath  region,  the  formulation  above 
onfy  models  the  potential  drop  portion  of  the  sheath,  and  does  not  attempt  to 
produce  the  potential  rise  of  the  inner  sheath.  In  addition,  Ecker’s  current 
constrictions  are  of  a  "macroscopic"  nature,  whereas  those  of  Reference  14  and  this 
work  are  "microscopic"  [Ref.  33]. 

Data  for  a  6,000" K  Nitrogen  plasma  were  used  to  test  the  equation  set  [Ref. 
21].  Producing  a  proper  solution  required  adjusting  the  initial  conditions  to  force  the 
curve  shapes  discussed  above.  Using  Equations  (2),  (3),  and  (4),  the  mean  free  path, 
Debye  length,  and  approximate  sheath  extent  are  calculated  as  X  »  9 . 352x10'^  m, 
Xp  »  1.690x10'^  m,  and  Xg  »  2.352x10*^  m  (this  assumes  a  drop  voltage  ot  100 
volts). 


20 


3.  >^nniiB«te  Fommlatton 


Reference  21  explores  the  above  equation  set  by  taking  advantage  of  the 
symmetry  among  the  equations,  and  introducing  two  parameters,  K*  and  K',  shown 
below. 


K*  s 

-K*  s 


It  can  be  shown  that  the  resulting  equations  can  be  manipulated  to  yield  a  single, 
ordinary  differential  equation  for  the  K’s  in  terms  of  the  electric  field.  The  resulting 
equation  can  be  scrutinized  for  two  distinct  temperature  regimes.  Note  that  while 
the  total  current  density,  J,  is  constant  in  a  steady,  one-dimensional  case,  the  K’s  can 
vary  and  will  in  turn  also  depend  on  the  degree  of  reactivity  of  the  plasma  (ft,),  i.e., 


S 

dy 


(19) 


Because  ion  diffusion  is  much  slower  than  electron  diffusion,  it  can  be  shown  that  the 


K’s  are  related 


K*  «  -K  —  (20) 

eD 

As  will  be  evident,  at  the  electrode  surface  the  K’s  are  equal  to  each  other  and  at  the 
undisturbed  plasma,  K~  s  o.  The  total  current  density  may  be  evaluated  from 
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J  -  cn.v 


(21) 


where  is  the  electron  drift  velocity  beyond  the  ambipolar  region  which  is  strialy 
a  hmction  of  E./N,  (i.e.,  of  the  ratio  of  undisturbed  electric  field  to  the  total  number 
density). 

a.  Effects  of  Tempemture  on  Anode  Constriction 

It  is  useful  to  investigate  the  overall  effects  of  temperature.  Since 
temperature  will  be  considered  constant,  it  comes  in  as  a  parameter  in  this 
formulation  whereas  charge  density  and  electric  field  remain  as  variables.  Intuitive 
arguments  will  be  introduced  which  suggest  that  the  electron  and  ion/neutral 
temperatures  play  a  rather  singular  role  in  determining  the  intrinsic  dimensionality 
of  the  problem,  (i.e.,  there  are  cases  when  the  geometry  of  the  current  lines  is  not 
necessarily  impressed  by  the  electrode  geometry).  Since  the  problem  is  described  by 
moderate  pressure,  largely  coUisional  sheaths,  the  ion  and  neutral  temperatures  are 
anticipated  to  remain  reasonably  equal.  Depending  on  the  gas,  the  electron 
temperature,  on  the  other  hand,  can  be  elevated  from  the  gas  temperature  at  the 
anode  where  actual  magnitudes  depend  on  the  local  value  of  E/N.  In  order  to  get 
a  perspective  on  the  effects  of  temperature,  we  shall  consider  two  extremes,  namely, 
the  case  where  the  electron  and  ion  temperature  are  the  same  (the  equilibrium  case) 
and  the  case  where  the  elearon  temperature  is  substantially  elevated  from  that  of 
the  ions/neutrals  (the  two-temperature  case). 
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(1)  Case  1:  T,  •  T.  •  T,  (Equmriurn) 


The  charge  densiti^  can  be  eliminated  by  combining  equations 
(5H9),  (17)  and  (18).  The  resulting  equation  can  be  shown  to  be 


2J 

eD 


kT 


•h  Li-  EE"-(Ey-I(-®-)TE" 
*  J 4'kT/ 


(22) 


If  the  electric  field  decreases  monotonically  from  the  wall  to  the  undisturbed  plasma 
(i.e.,  from  Eo  -*  E,),  then  as  y  -*  «,  E  -►  E«  E’  -»  0,  E”  -»  0. 

So  that  in  equation  (22)  above  the  "outer  solution"  becomes: 


K*  -  _  (23) 

eD. 

Now  this  represents  an  acceptable  solution  from  a  physical  point  of  view.  Moreover, 
asy-*  «  , 

K  •  D.(KT  *  0  (24) 

whidi  is  also  acceptable  for  an  equilibrium  situation  at  the  undisturbed  plasma. 
Results  [Ref.  21]  are  shown  in  Figure  5  for  the  case  of  nitrogen  at  6000”K  using  an 
approximate  electric  field  distribution. 
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Figure  5  -  Electric  Held  and  Species  as  a  Function  of 
Approximation  Using  a  Shaped  Electric  Field  and  Iso 


2* 


(2)  Case  IJ:  T,  >>  T,  »  T,  (Two'Tenq>erature) 

In  this  case  the  same  procedure  as  before  yields  the  following 
equation  where  terms  divided  by  T,  have  been  drqf^d  udien  conqwred  to  their 
counterparts  divided  T,. 


^  c'  _  2EJ 
^  ^  kTo: 


EE"  ♦  — E'E'  -  ilE'  (25) 
kT,  eE  e 

Assuming  the  same  monotonic  decrease  as  before  for  the  electric  field  from  the  wall 

to  the  plasma  proper,  as  y  -»  «,  E  -  E„  E’  -►  0,  E”  0. 

Then  the  outer  solution  becomes 


dK- 

dy 


2eEJ 

kT,eD. 


with  h^>0 


(26) 


Or,  K*  (constant)  y  constant,  and  keeps  increasing  with  y. 

This  is  not  the  prc^r  outer  solution  for  the  one-dimensional,  equilibrium 
plasma  that  we  seek  because  the  net  ionization  rate  continues  to  increase  well  inside 
the  plasma  prc^r  where  conditions  should  saturate,  yielding  a  constant  electric  field. 
Therefore,  as  formulated.  Case  n  is  not  amenable  to  a  one-dimensional  solution. 
References  14  and  21  show  how  this  case  can  be  analyzed  under  a  multidimensional 
£q>proach.  These  references  also  discuss  a  method  for  describing  the  electron 
tenq)erature  as  a  function  of  E/N,  then  how  to  couple  a  simplified  energy  relation 
ti^iidt  satisfactorily  describes  a  two-temperature  plasma.  The  necessary  ingredient 
to  make  equation  (26)  approach  zero  beyond  the  decrease  of  E  to  E.  is  to  allow  J 
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to  fan  out  as  indicated  in  Figure  6.  Thus,  in  equation  (26),  the  product  "EJ"  can 
bring  down  the  charge  production  rate  to  arbitrarily  low  values.  Alternatively,  it  is 
possible  to  explore  techniques  of  bringing  the  electron  temperature  down  to  be  in 
closer  equilibrium  with  the  ions  and  neutrals.  Transpiration  cooling  is  one  such 


means. 


Figure  6  -  Two-Dimensional  Model  of  Current  Paths  Showing  Periodic  Structure. 
Thermal  Instabilities  and  Inhomogeneities  Would  Favor  One  Site  Over  Others  and 
a  Single  Macroscopic  Constnction  May  Then  Be  Produced.  [Refs.  14,  34 J 
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b.  SMmbif  to  Vacuun  Ak  Phenomena 

IiotaMUty  phencnnena  observed  in  vacuum  arcs  [Ref.  35]  are  veiy 
similar  to  those  d»eived  in  self-field  thrusters  [Ref.  12].  After  the  establishment  of 
the  current,  the  anode  region  operates  in  a  v^r  that  issues  from  the  electrodes. 
In  vacuum  arcs.  Miller  diaracterizes  the  anode  region  as  operating  in  one  of  five 
distinct  modes,  ranging  frmn  a  passive,  low  current  mode  to  a  high  current,  fully 
devetoped  spot  mode  [Ref.  36].  Given  the  simiUrities  mentioned  above,  vacuum  arc 
anode  research  should  be  helpful  in  the  understanding  of  MPD  thruster  transition 
to  the  anode  spot  mode.  Existence  diagrams  after  Miller  [Ref.  36]  are  shown  in 
Figure  7,  which  divide  operating  modes  into  regions  as  a  function  of  anode  current 
versus  electrode  geometry.  Figure  7  shows  the  transition  from  glow  to  spot  mode. 

Anode  ^t  formation  at  high  currents  is  clearly  a  factor  in  limiting  anode 
lifetime.  Various  pheimmena  have  been  related  to  anode  spotting.  Hu^l  [Ref.  12] 
relates  the  transition  to  spotting  mode  to  an  increase  in  J^/ih  above  a  critical  level. 
A  separate  factor  connected  with  the  spot  mode  is  surface  temperature  of  the  anode. 
Rich,  etal.,  [Ref.  37]  show  that  anode  spotting  is  preceded  by  a  luminous  "footpoint" 
and  followed  by  local  melting  prior  to  spot  formation.  Separately,  Schuocker  [Ref. 
38]  finds  a  connection  between  spotting  initiation  and  the  factors  of  anode 
evaporation  and  magnetic  constriction  in  vacuum  arcs  with  high  currents. 
Ejq)erimental  investigations  must  be  performed  to  see  if  the  above-mentioned 
vacuum  arc  criteria  apply  to  self-field  thrusters. 
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ARC  CMRCMT  (M> 


Fi^e  7  •  Anode  Discharge  Modes  as  a  Function  of  Current  and  Gap  Length.  [Ref. 


C  COMPUTER  CODE 


^ 


Ratter  than  using  linear  approximations  to  equations  (12H16),  the  nonlinear 
set  was  used,  with  initial  conditions  adjusted  in  an  attempt  to  produce  observed 
electric  fields  from  probe  data.  First  and  second  space  derivatives  of  the  electric  field 
were  used  as  diagnostic  checks  to  ensure  computed  output  was  reasonable.  Initial 
conditions  computed  from  the  approximate  formulae  in  Reference  21  were  used.  The 
equation  set  above  presents  a  difficult  problem  for  two  reasons,  nonlinearity  and 
multiple  time  constants.  The  species  number  density  equations,  (12)  and  (13),  both 
contain  a  nonlinear  term,  each  with  a  time  constant  of  its  own.  In  addition,  the 
electric  field  equation,  (14),  adds  a  possible  third  time  constant.  This  constitutes  a 
"stiff*  set  of  equations.  Attempts  were  made  to  solve  the  set  with  the  data  discussed 
above,  using  Gear’s  method  of  backward  differentiation,  in  hopes  that  the  variables 
would  change  slowly  enough  with  each  iteration  to  render  a  convergent  iterative 
process.  As  described  in  Reference  39,  if  some  reactions  are  slow  and  others  fast 
among  a  set  of  coupled  equations,  the  fast  ones  will  control  the  stability  of  the 
method.  This  is  addressed  in  the  DGEAR  program  available  from  the  International 
Mathematical  &  Statistical  Library  (IMSL).  The  latter  software  contains  an  Adams 
predictor-corrector  method,  as  well  as  Gear’s  method,  which  is  well  known  for  its 
success  at  solving  stiff  equation  sets.  The  DGEAR  software  allows  for  a  choice  of 
fimctional  or  chord  iteration  methods,  as  well  as  a  choice  of  Jacobian  matrices.  A 
more  detailed  discussion  of  this  software  can  be  found  in  Reference  39  and  in  the 
IMSL  Ubraiy.  [Ref.  39] 
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D.  COMPUTATIONAL  RESULTS 


Numerous  computer  runs  were  completed  using  the  initial  conditions  taken 
from  Reference  21.  In  addition,  data  for  the  ionization  coefficient  v ,  Figure  8,  was 
taken  from  References  40  and  41. 


Figure  8.  Ionization  Coefficient  v  as  a  Function  of  E/N  for  Nitrogen  for  Various 
Vibrational  Temperatures  (Refs.  40,41). 

Various  combinations  of  initial  conditions  and  ionization  coefficient  were  used.  As 
mentioned  above,  the  electric  field  and  its  space  derivatives  were  used  as 
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(tiagpostk/reasonability  cfaedcs  on  the  outinit.  Individual  as  well  as  multiple 
ccHiqpiiter  nins  were  attempted  to  model  the  sheath  region.  Nonlinearities  in  the 
equation  set  are  dearly  seen  in  Figure  9.  The  ion  number  density  does  not  reach  that 
of  the  electrons,  and  the  latter  peculation  growth  rate  continues  to  grow  without 
bound.  The  shape  of  the  electron  population  curve  is  very  sensitive  to  its  initial  value. 
As  shown  in  Figure  9,  the  latter  population  has  too  hi^  a  growth  rate  when 
compared  to  the  ion  peculation,  and  the  latter  does  not  "catch  up".  Increasing  the 
initial  value  of  fl,  flattens  out  this  curve  to  a  reasonable  shape.  Above  an  initial 
value  of  0.06,  he>wever,  the  plot  of  "dips"  after  a  certain  distance 

and  then  continues  to  increase  as  expected.  This  gives  an  sqcfodmate  ujc^r  value 
for  this  initial  value.  To  avoid  instabilities  like  this,  small  "slices"  were  taken  of  the 
ouqmt  after  a  small  number  of  integration  steps  and  multiple  runs  were  used  to  form 
a  "cut  and  paste"  plot  of  the  region.  When  a  reasonable  plot  shape  was  produced,  the 
value  of  iemization  coeffident  was  varied  in  the  "slices"  to  attempt  to  produce  the 
required  end  values  for  electric  field  and  spedes  population.  Both  multiple  and 
equilibrium  values  for  the  ionization  coeffident  were  used.  When  the  data  showed 
sips  of  instability  and  failure  to  follow  the  required  forms  of  Fipre  4,  a  "slice"  was 
made  in  the  data  stream,  and  the  data  points  from  this  point  used  to  start  a  new 
computer  run.  This  approach  was  taken  in  the  hope  of  avoiding  singularities  in  the 
integration  from  anode  surface  to  ambipolar  region.  In  addition  to  the  diapostic 
checks  shown  in  Fipre  4,  an  additional  data  check  is  provided  by  the  transition  from 
the  sheath  to  the  ambipolar  region. 
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Figure  9.  Species  Number  Density  plots  for  Individual  Computer  Run,  Showing 
Divergent  Tracks  for  Ion  and  Electron  Populations,  and  Effect  of  Nonlinearity. 


As  shown  in  Figure  5,  the  species  number  densities  are  equivalent  in  this  region,  as 
are  their  change  rate.  Thus,  setting  Equations  (12)  and  (IS)  equal  to  each  other  and 
solving  for  f1  yields  a  value  of  0.5  in  the  ambipolar  region.  As  indicated  in  Figures 
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10>1U  the  output  produces  the  desired  plot  slopes  for  electric  field  and  species 
number  density.  However,  the  number  density  plots  cross  long  before  approaching 
the  required  value  of  OJ.  In  addition,  neither  electric  field  nor  species  number 
density  approaches  an  equilibrium  value  or  shows  sign  of  levelling  off.  Apparently, 
the  multiple  time  constants  and  nonlinear  portions  of  the  number  density  equations 
combine  to  create  a  seemingly  intractable  system.  Solutions  for  this  system  may  be 
possible  for  specific,  individual  initial  condition  sets,  but  the  problem  does  not  appear 
amenable  to  this  approach  in  general.  A  one-dimensional  system  such  as  this  may  be 
better  described  through  the  approach  of  boundary  layer  theory  or  nonlinear 
dynamics  and  chaos.  Given  the  effort  and  difficulty  involved  in  the  latter,  a  one¬ 
dimensional  approach  such  as  that  modelled  above  does  not  appear  useful.  A 
combination  of  one-  and  two-dimensional  modelling  would  appear  to  be  more  useful, 
as  discussed  in  Reference  14.  A  one-dimensional  model  may  be  useful,  but  only  in 
an  approximation  approach,  with  a  shaped  electric  field,  such  as  that  used  in 
Reference  21. 
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Figure  10.  Electric  Field  as  a  Function  of  9,  Distance  From  Anode,  Using 
Equations  (12)*(16). 
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ionization  (with  some  help  from  the  tail  of  the  Maxwellian  distribution  of  electron 
energies),  what  results  is  a  breakdown  voltage  aj^reciably  below  the  ionization 
potential.  This  then  could  be  an  explanation  for  the  low  voltage  breakdown 
observations  [Ref.  42].  Qearly,  gases  with  low  ionization  potentials  and  lots  of 
atomic  electron  energy  levels  are  preferred  (such  as  cesium  and  barium)  but  low- 
voltage  breakdown  has  been  observed  with  most  gases. 

The  increase  of  the  anode  fall  voltage  above  the  ionization  potential  has  been 
related  to  the  electron  Hall  parameter,  since  a  reduction  of  this  parameter  decreases 
that  voltage  and  corresponding  losses.  Control  of  the  local  magnetic  field  through 
the  use  of  and  array  of  permanent  magnets  as  well  as  implementation  of 
transpiration  cooling  (which  increases  the  electron  collision  frequency)  have  both 
yielded  some  encouraging  results.  Because  the  anode  fall  also  scales  up  with  jV^, 
it  is  conceivable  that  current  inhomogeneities  and  plasma  instabilities  which  are 
reflected  in  this  parameter  are  in  the  picture  as  well.  [Ref.  44] 

In  summary,  any  possible  reduction  of  the  anode  fall  voltage  will  hinge  on  a 
thorough  understanding  of  the  anode  region,  with  its  associated  sheath  and  ambipolar 
regions,  where  electron  temperature  effects,  ionization  effects,  and  magnetic  field 
effects  play  a  pivotal  role.  If  transpiration  cooling  is  present,  then  additional 
phenomena  of  fluid-dynamic  nature  may  come  into  play.  Experimental  observations 
with  atmospheric  discharges  indicate  the  possible  presence  of  convective  effects  at 
the  anode.  [Ref.  45] 
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IV.  TRANSPIRATION  COOUNG 


TnuK^nration  cocking  of  the  anode  has  often  been  promoted  as  an  attractive 
means  oi  recovering  a  large  portion  of  the  power  deposited  there.  Additionally,  the 
onset  of  melting  may  be  minimized  or  even  avoided  by  active  anode  cooling.  Ridi, 
etal.,  related  high  anode  ten^ratures  to  anode  spotting  [Ref.  37].  Similarly,  Park 
and  Choi  showed  that  low  thermal  difhision  leads  to  erosion  and,  consequently, 
anode  damage  [Ref.  46].  Active  anode  cooling  via  transpiration  is  one  means  of 
ensuring  high  thermal  diffusion  and  extending  anode  lifetime.  Early  work  by 
Schoeck,  et  al.,  [Ref.  47]  showed  that  up  to  80%  of  the  energy  deposited  in  the 
ano^  is  recoverable  via  transpiration  cooling.  While  this  study  used  non-convective, 
high-intensity  argon  arcs,  it  is  reasonable  to  assume  that  this  effect  would  apply  in 
MFD  arcs  uang  other  propellants.  Although  this  cooling  method  has  not  been 
studied  for  incorporation  in  plasma  thrusters  for  some  time,  it  has  been  recently 
considered  as  a  means  of  cooling  the  fuselage  of  the  National  Aerospace  Plane 
(NASP).  Plasma  thruster  designs  could  undoubtedly  gain  from  this  database,  and  due 
consideration  should  be  given  to  this  cooling  approach  for  the  anode. 

For  a  given  mass  transfer  flow  rate,  the  heat  flux  reduction  to  a  surface  is 
inversely  proportional  to  the  molecular  weight  of  the  injected  gas.  Use  of  the 
prqjellant  as  coolant  as  well  as  fuel  would  eliminate  the  need  for  additional  tankage 
and  pumps,  simplifying  the  design  considerably.  Lithium  has  been  considered  to  be 
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the  pn^Uant  of  choice,  primarily  because  of  its  low  molecular  weight,  its  favorable 

ionization  potential,  and  its  low>volume  tankage  properties.  It  has  a  relatively  low 

first  ionization  potential  of  5.4  eV  and  a  high  second  ionization  potential  of  7S.6  eV. 

This  single  ionization  potential  range  of  over  70  eV  compares  to  approximately  20 

eV  for  Cesium  and  27  eV  for  Potassium  [Ref.  48].  This  provides  a  broad  temperature 

range  within  which  only  single  ionization  will  occur.  Large  ten^ratures  must  be 

readied  within  the  gas  before  double  ionization  occurs.  As  the  gas  temperature  is 

increased  several  thousand  degrees  Kelvin,  it  undergoes  ionization  and  disassodation. 

Thermal  energy  deposited  can  be  recovered  through  nozzle  expansion  at  the  exit. 

However,  residence  time  of  the  gas  is  not  long  enough  to  ensure  recombination. 

Thus,  the  energy  invested  in  ionization  and  disassodation  will  be  lost.  [Ref.  49] 

lithium  has  been  shown  to  produce  specific  impulse  figures  in  excess  of  7,000 

seconds  at  70%  effiden^  in  steady-state  thrusters,  [Ref.  50]  whereas  all  other 

propellants  have  been  limited  to  less  than  3,000  seconds  spedfic  impulse  at  less  than 

40%  effiden<7  [Ref.  6].  Subramaniam  has  concluded  that: 

..regenerative  cooling  of  anodes  (at  the  specific  impulse  values  in  the  MPD 
regime)  is  possible  only  with  hydrt^en  or  with  alkali-metal  propellants,  notably 
lithium.  In  the  latter  case,  the  ideal  anode  operating  mode  would  be 
evaporation  and  ionization  of  the  propellant  on  the  porous  or  wetted  anode 
surface,  resulting  in  increased  ion  current  fraction,  reduced  anode  voltage  fall 
and  utilization  of  part  of  the  anode  loss  energy  [Ref.  51]. 

liquid  coolants,  as  well  as  redudng  storage  requirements,  offer  the  advantage  of 

providing  latent  heat  of  vaporization  for  energy  disposal.  However,  design  problems 

can  occur  if  the  liquid  is  allowed  to  vaporize  within  the  porous  structure.  Problems 
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arise  due  to  the  abrupt  increase  in  pressure  gradient  as  the  coolant  vi^rizes.  Since 
coolant  flow  generalfy  have  three>dimensional  characteristics,  the  flow  will  be 
diverted  around  the  vapor  bubbles  and  hot  ^>ots  often  develop.  The  technical 
practicality  of  using  molten  lithium  to  cool  a  porous  tungsten  anode  would  seem  to 
be  beyond  current  technology.  On  the  other  hand,  the  products  of  decomposition  of 
hydrazine  (gaseous  hydrogen  and  ammonia)  have  proven  to  be  efflcient  and  practical 
coolants  [Ref.  S2]. 

Given  the  performance  figures  above,  using  an  auxiliary  coolant  gas  even  with 
high  molecular  weight  (e.g.,  NH„  Nj,  CH4,  etc.)  which  could  serve  as  a  propellant 
once  released  firom  a  porous  hot  tungsten  anode  surface  would  seem  to  more 
practical,  vice  dealing  with  molten  lithium.  Experimental  studies  would  be  needed 
to  compare  the  approaches.  Kuriki  and  Suzuki  performed  erqreriments  with  a  quasi¬ 
steady  MPD  thruster  to  study  the  effect  of  anode  gas  injection  (Argon).  At  high 
currents  of  up  to  10  kA,  increases  in  thrust,  specific  impulse,  and  flow  discharge 
stability  were  observed  [Ref.  53]. 

There  is  some  question  as  to  the  likelihood  of  current  constriction  resulting 
from  anode  gas  injection.  In  such  a  case,  swirling  or  circulating  the  propellant  gas 
would  help  to  move  any  footpoints  that  developed  around  the  anode  surface  and 
prevent  them  fi'om  becoming  fully-developed  spots.  Additionally,  an  applied 
magnetic  field  could  serve  to  circulate  the  footpoints  as  well.  The  unique  advantage 
of  transpiration  cooling  hinges  on  providing  effective  anode  cooling  while  supplying 
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liot  im^llant,  but  the  real  ben^t  will  depend  on  how  small  the  amount  of  a>olant 
letpiired  really  will  be. 

Trat^inration  cooling  has  proven  to  be  as  desiraUe  as  it  is  challenging.  It  is 
ctmqriicated  to  implement,  with  associated  reliability  im>blems  and  difficulty  of 
analytical  predictiom.  While  the  productkm  of  thidcer  boundary  la^rs  is  largely 
ineffective  against  the  electron  flux  heating,  the  cooling  itself  is  nrnst  effidem  and  a 
substantial  fraction  of  the  energy  transferred  to  the  anode  is  recoverable.  The 
arguments  of  Chapter  Three  indicate  that  a  reduction  of  the  electron  temperature 
in  the  anode  would  have  the  desirable  effect  of  reducing  the  initial  current  spotting 
whidi  can  be  conjectured  to  be  the  path  that  leads  to  anode  arc  ^ts.  This  electron 
temperature  reduction  can  be  deme  m(»t  effectively  by  polyatomic  pses  (which  have 
a  higb  5-loss  factor)  emanating  from  the  anode  surface  [Ret  54]. 

The  arguments  relating  to  tranqpiratimi  cooling  might  be  summarized  as 
ftdlows: 

Favorable  Outcomes 

•  No  separate  cotding  mechanism  for  anode  required, 

•  Adds  "hot"  prqiellant  to  exhaust  "recovering”  most  of  the  electrical  power  loss  to 
the  anode, 

•  Quenches  T.  thus  likely  to  postpone  anode  spotting  ^  reduce  the  heating 
associated  with  the  electron  thermal  energy  (SkT^Ze), 

•  Reduces  bulk  convective  heating, 

•  Rechices  the  local  electron  Hall  parameter  by  increasing  the  collision  frequency. 
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Altows  fcff  some  rediatkm  cooling  from  the  hot  tungsten  surface  (abmit  12K) 


watts/cm*  at  TSOBtK  (Ret  51]), 

•  H^dn^en/ammonia  gases  flowing  through  hot  porous  (sintered)  tungsten 
rqwesent  a  cmiqtatitde,  proven  technology. 
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•  Mi^  decrease  the  electrical  conductivity  in  the  anode  region, 

•  May  destabilize  the  icmization  processes  in  the  sheath  and  bring  about  significant 
fluctuati<His  in  the  current, 

•  IKsrupts  "cathode  jet"  in  front  of  the  anode  with  unpredictable  consequences, 

•  Introduces  propellant  udiidi  may  not  be  hot  enough,  not  icmized  enmigh,  or  not 
in  the  proper  place  fw  J  xS  acceleration, 

•  Transpiratton  cooling  throu^  a  porous  (tungsten)  anode  is  a  difficult  des^ 
problem. 


V.  CONCLUSIONS  AND  RECOMMENDATIONS 


nasma  thrusters  offer  distiiict  advantages  in  terms  of  payload  delivered  for 
interidaiietaiy  missi<Kis,  as  well  as  for  mlntal  transfer.  A  recent  comparison 
completed  by  Chmieiri,  Kelfy,  and  Jahn  shows  a  mass  savings  of  65  tons  for  an 
mbital  tramfer  from  low  Earth  orbit  to  geosynchronous  Earth  orbit  using  a 
quasisteatfy  MFD  thruster  as  q^)Osed  to  an  advanced  chemical  thruster.*  This 
siq>erior  performance  comes  at  the  e]q)ense  of  low  thrust-to-weight  ratio  and  long 
transit  time.  However,  given  the  large  cargo/logistic  requirements  of  a  manned 
interplanetary  mission,  delivery  of  payload  must  be  maximized.  Thus,  further  work 
to  diaracterize  more  fully  thruster  behavior  and  anode  contributions  in  particular  are 
certainly  warranted.  [Ref.  55] 

The  "cut-and'paste"  method  used  to  generate  Figures  10  and  11  is  not  of 
practical  use  as  a  modelling  method,  due  to  the  large  effort  involved.  It  did  produce 
the  expected  electric  field  and  species  number  and  current  density  plots  near  the 
anode,  but  failed  to  produce  the  entire  sheath  out  to  the  ambipolar  region.  The 
nonlinearity  of  the  equation  set  led  to  a  quickly  deteriorating  solution.  A  more 
practical  approach  using  nonlinear  dynamics  and/or  chaos  must  be  developed  to 
model  the  sheath  numerically. 


"This  assumes  a  specific  impulse  of  2,000  seconds,  600  kW  of  input  power,  and 
a  270-day  transit  time. 
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Depending  on  the  prt^Uant  mass  frt^on  used  for  cooling,  the  transpiration 
scheme  discussed  above  presents  some  rather  unique  advantages.  A  hot  anode 
whidi  uses  only  a  miall  amount  of  prt^Uant  for  cooling  need  not  be  penalized  for 
any  lost  thrust  If  in  addition,  we  increase  anode  lifetime  by  delaying  the  formation 
of  anode  arc  ^ts,  then  the  scheme  is  all  the  more  desirable.  A  decrease  of  the 
electron  temperature  in  the  vicinity  of  the  anode  may  bring  about  a  more 
homogeneous  flow  of  current  and  a  reduction  in  the  heating  effect  associated  with 
the  high  electron  kinetic  energy.  Recovery  of  the  heat  deposited  at  the  anode  would 
be  most  important  if  the  propellant  fraction  is  high.  In  such  case,  nozzle  expansion 
of  the  hot-propellant/coolant-gas  might  be  implemented. 

Means  of  limiting  anode  losses  through  decreasing  anode  fall  voltages  were 
discussed,  including  the  control  of  the  local  Hall  parameter  and  the  implementation 
of  thermionic  arc  breakdown.  The  electrical  conductivity  (of  a  nonreacting  plasma) 
could  possibly  decrease  as  a  result  of  transpiration  cooling  and  this  might  increase 
the  anode  fall  voltage. 

Additional  work  needs  to  be  done  in  the  following  areas: 

•  Investigate  effectiveness  of  nonlinear  dynamics  and  chaos  in  solving  sheath 
equation  set, 

•  Incorporate  adequate  one-  or  two-dimensional  sheath  modelling  in  quasisteady 
MPD  numerical  codes, 
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•  Investigite  the  role  that  fluid  dynamic  effects  play  in  MPD  thruster  anode 
disdiarges, 

•  Investigate  the  effect  of  tran^iration  cooling  on  current  and  plasma  stability,  as 
well  as  on  thruster  performance  and  lifetime, 

•  Determine  effectiveness  of  transpiration  cooling’s  increase  of  the  collision 
frequent  parameter, 

•  Compare  performance  of  gaseous  propellant/coolants  versus  hybrid  designs  with 
lithium  propellant/gaseous  coolant, 

•  Determine  if  required  percentage  of  propellant  gas  as  coolant  is  practical  (e.g.,  less 
than  10%), 

•  Investigate  effect  of  surface  imperfections  as  focal  points  for  current  constrictions 
and  as  precursors  to  anode  spotting. 
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The  following  software  includes  the  calling  program,  SHBATH.  its  two 
subroutines,  FCSU  and  TPOT,  and  the  DGBAR  integrator.  The  latter  is  quite 
extensive  in  length  and  includes  ten  subroutines,  including  the  following:  DORST, 
DGRCS,  DORPS,  DGRIN,  LOCATF,  lOBLMP,  UtQTlB,  DBRTST,  DOBTIO,  and  USPKD.  A 
detailed  discussion  may  ba  found  in  lUSL  literature  or  Reference  39. 


Program  Sheath  OOOOIO 

C  000020 

C . Calling  program  for  DGKAR  integrator.  Initial  conditi^is  are  000030 

C  input  via  RSAD  statesients  and  keyboard  entry.  Output  is  to  data  000040 

C  files  via  the  DORST  subroutine.  Diagnostic  check  of  output  via  000050 

C  Figure  4  printed  to  data  file  from  DORST  subroutine.  Consult  000060 

C  DOBAR  caanients  for  variable  descrlpticns  not  listed  below.  000061 

C .  000062 

RBAL  B,K,BPS,TI,BF,BFINP,DZ,DB,innHF,Cl,Kl,A  000070 

INTEGER  N.MBTH,MITBR, INDEX, imCd) ,IER, STEP  000080 

RBAL*8  X,H,Y(5) ,XBND,TOL,WK  000090 

EXTERNAL  YDOT,  FCNJ,  DOBAR  000100 

CGMMON/CONST/B , K, BPS , TI , BF , BFINF, NINF , DI . DB , VINF ,C1,K1,A  000110 

C  000120 

C . Constants  000121 

C  Cl  and  Kl  are  constants  describing  the  ionisation  coefficient.  000122 

C  They  are  taken  from  the  data  plotted  in  Figure  8.  The  000123 

C  coefficient  is  equal  to  the  nondimensionalised  electric  potential  000124 

C  raised  to  the  Kl  potier  and  then  multiplied  by  Cl.  000125 

C  In  this  way,  the  ionization  coefficient  is  allowed  to  vary  in  000126 

C  proportion  to  the  strength  of  the  electric  potential.  000127 

C .  000128 

WRITE (*,*)' Input  value  for  Cl  (format  6E3):'  000129 

READ(*,*)C1  000130 

WRITE (*,*) Cl  000131 

WRITE (*,*)' Input  value  for  Kl  (format  6B3):'  000132 

RBAD(*,*)K1  000133 

WRITE (*,*)K1  000134 

Initial  conditions  for  species  number  density,  electric  potential  000135 
euid  species  current  density  are  now  input  (ni,ne,E, je, ji) .  000136 

.  000137 

WRITE (*,*) 'Input  values  for  y(l)  through  y(5)  (format  5(6E3}):'  000138 

READ(*,*)y(l)  ,y(2)  ,y(3)  ,y(4)  ,y(5)  000139 

WRITE (*,*)y(l) ,y(2) ,y(3) ,y(4) ,y(5)  000140 

Following  constauits  are  for  plasma  described  in  Reference  21  000141 

(6,000  K,  Init  E»20,000  V/m,  Final  E-1,200  V/m)  000142 

B-1.6B-19  000143 

K«1.38B-23  000150 

BPS>8.854E-12  000160 

TI«6B3  000170 

BF«2B5  000180 

BFINF«1.2E4  000190 

DI«1.724B-4  000200 

DB-1.724E-1  000210 

VIO  -  2.B6  000215 

VINF  >  4.93B-7  000220 

C  000230 
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onoo  ooo  0^0 


C 


A  is  plj 


cAisrsctsristic  Isngth  whidti  shows  potsatisl  drop. 


000240 


C  A«  ( (BPS*1P)  /  (I^HIinP) )  .1.1071-6 

A  .  1.107B-S  -  000250 

X  -  0.01  000270 

XWP  a  10.  000280 

H  ■  lS-6  000290 

TOli  a  11-6  000300 

MBTH  a  2  000305 

MXTBR  a  1  000310 

I8DBX  a  1  000320 

H.5  000330 

imcd)  a  5  000340 

me  a  18000.  000350 

ISR  a  0  000360 

OPIN(aNlTa8.FILBa'SHBATH.DAT'  ,STATaSa'iniiaK>IIN' )  000370 

CAU.  DOBAR2(N.YDOT,PCNJ,X.R«Y.XBHD,TOL.MBTO, MXTBR,  INDBX,  IKK,  WK,  000380 
-t-IBR.STBP)  000390 

DO  3  Ia0,N  000400 

DO  2  JaO.lOO  000410 

MRrrB(*.*)J,Y(I)  000420 

NRITB(8,1) J.Y(I)  000430 

1  FORMAT(T2,F5.1,5(5X,D9.2) )  000440 

2  COMTIRDB  000450 

3  COrriMOB  000460 

1IR1TB(«,*) 'Total  Stops  a  ' ,STBP, 'Final  Step  Size  «  ' ,H,  000470 

Error  Code  .  ',IBR  000480 

CXiOSB(T]llITa8)  000490 

STOP  000500 

BMD  000510 


DOMIY  SUBRODTINB  FCNJ 


SUBROOTIMB  FCNJ(N,X,Y,PD)  1 
IIITB6BR  N  2 
REAL  Y(N}  ,PD(M,N)  ,X  3 
RETURN  4 
END  5 


SUBROOTINB  YDOT 


SUBROUTINE  YDOT(N,X, Y, YPRIMB,eprinie,epriine2) 

RBAL*8  X,  Y(5)  ,YPRIHB(5)  ,RUI,epriaie.epriilie2 
REAL  B,K,BPS,TI,BF,BFIRP,NINF,DI,DE,VINF,C1,K1,A,B1,B2,B3,B4 
COHMON/CemST/E ,  K,  BPS ,  TI ,  BF ,  BFINF ,  NINF,  DI ,  DB ,  VIO ,  VINF ,  Cl ,  K1 ,  A 
VI  a  Cl  *  (y(3)**Kl) 

VIT  a  VI  /  VIO 

Following  constants  are  the  bracketed  values  in  Equations  12-16. 
A  is  left  as  a  variable. 


B1  a  {(B*BPS)/(K*TI))  *A 
Bla  3.86B5  *  A 
C  B2  a  { (B*EFINF) / (K*TI) )  *  A 
B2  a  2 . 32B4  *  A 

C  B3  a  ( '(B*NINF)  /  (EF*BPS) )  *  A 
B3  a  9.04B5  *  A 

C  B4  a  ((VINF*K*TI)/(E*DB*BFINF))  *A 

B4  a  2.62B-21  *  A 

C  Alpha  a  2 -body  recombination  coefficient  (fm.  Laser  Kinetics 
C  Handbook  (AFNL-TR-74-216,  1974))  (cm3/sec) 

Alj^ha  a  9.e-8 
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C . 

C  PIVI  FIRST  OSDIR  IQOKnCMS  -  Iquatloos  12-16 


Hi 

TFRIMKl)  -  (B  *  Yd)  *  Y(3))  -  Y(5) 

M* 

YFRI]IB(2)  •  -(B  *  Y(2)  *  Y(3))  -t-  Y(4) 

B 

YFRI1IB(3)  •  B3  *  (Y(l)  •  Y(2)) 

YPRI1IB(4)  >  -B4  *  Y(2)  «  (VIT  -  (ALPHR  *  Y(l))) 

ji 

YPRIMB(5)  -  B4  *  Y(2)  *  (VIT  -  (ALPHA  *  Y(l))) 

Diagnostic  Checlc  of  first, ssccad  derivativss . 

aRrina  «  y(l)  -  y(2) 
aprinr2  >  yprimsd)  -  ypriins(2) 

C 

RBTDRN 

BHD 
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c 

IMSL  ROOTISB  NAMB 

-  DGEAR 

DGEAOOIO 

c 

DOBA0020 

c- 

•modified  to 

return 

#  of  steps  via  variable  "step*  in  siibroutine  call  * 

c 

DOBA0040 

c 

COUPDTBR 

-  IBM/DOUBLE 

D(»A0050 

c 

DOBA0060 

c 

LATEST  REVISION 

•  NOVEMBER  1,  1984 

DGEA0070 

c 

DGEA0080 

c 

PORPOSB 

-  DIFFERENTIAL  E(^TION  SOLVER  •  VARIABLE  ORDER  DGEA0090 

c 

ADAMS  PREDICTOR  CORRECTOR  METHOD  OR 

DCTAOlOO 

c 

(»ARS  METHOD 

DGBAOllO 

c 

DGBA0120 

c 

USAGE 

-  CALL  DGEAR  (N,FCN,FCMJ,X,K,Y,XEND,TOL,MBTH, 

DGBA0130 

c 

MITER, INDEX, INK, NK, lER) 

DaEA0140 

c 

DGBA0150 

c 

ARGUMENTS 

N 

•  INPUT  NUISER  OF  FIRST-ORDER  DIFFERENTIAL 

DGBA0160 

c 

EQUATIONS . 

DOBA0170 

c 

FCN 

-  NAME  OF  SUBROUTINE  FOR  EVALUATING  FUNCTIONS. 

DGEA0180 

c 

(INPUT) 

DOBA0190 

c 

THE  SUBROUTINE  ITSELF  MOST  ALSO  BE  PROVIDED 

DOBA0200 

c 

BY  THE  USER  AND  IT  SHOULD  BE  OF  THE 

DGBA0210 

c 

FOLLOWING  FOIOf 

DGBA0220 

c 

SUBROUTINB  FCN  (N,X, Y, YPRIME) 

DGEA0230 

c 

REAL  X,  Y  (N)  ,  YPRIME  (N) 

D6BA0240 

c 

DOBA0250 

c 

DQBA0260 

c 

DaBA0270 

c 

FCN  SHOULD  EVALDATE  YPRIME  (1)  ,  . . . ,  YPRIME  (N) 

DGBA0280 

c 

GIVEN  N,X,  AND  Y(l) ,  . . . ,  Y(N)  .  YPRIME  (I) 

DGBA0290 

c 

IS  THE  FIRST  DERIVATIVE  OF  Y(I)  WITH 

DGBA0300 

c 

RESPECT  TO  X. 

D(XBA0310 

c 

FCN  MUST  APPEAR  IN  AN  EXTERNAL  STATEMENT  IN  DOEA0320 

c 

THE  CALLINB  PROGRAM  AMD  N,X,Y(1)  ,  .  .  .,Y(N) 

DaBA0330 

c 

MUST  MOT  BE  ALTERED  BY  FCN. 

DOBA0340 

c 

PCNJ 

•  NAME  OF  THE  SUBROUTINE  FOR  COMPUTING  THE 

DaBA0350 

c 

JAC(»IAN  MATRIX  OF  PARTIAL  DERIVATIVES. 

DaEA0360 

c 

(INPUT) 

DOBA0370 

c 

THE  SUBROUTINE  ITSELF  MUST  ALSO  BE  PROVIDED 

DGEA0380 

c 

BY  THE  USER. 

DCaCA0390 

c 

IF  MITBR-1  IT  SHOULD  BE  OF  THE  FOLLOWING 

DGBA0400 

c 

FORM 

DGEA0410 

c 

SUBROUTINE  FCNJ  (N,X,Y,PD) 

DOBA0420 

c 

REAL  X,Y(N)  ,PD(N,N) 

DGEA0430 

c 

DGBA0440 

c 

DGEA0450 

n 

FCNJ  MUST  EVALUATE  PD(I,J)  ,  THE  PARTIAL 

DOBA0460 

c 

DERIVATIVE  OF  YPRIME  (I)  WITH  RESPECT  TO 

DGEA0470 

c 

Y(J)  ,  FOR  I«1,N  AND  J-1,N. 

DGBA0480 

c 

IF  MITER-  -1  IT  SHOULD  BE  OF  THE  FOLLOWING 

DGEA0490 

c 

FORM 

DGEA0500 

c 

SUBROUTINE  FCNJ  (N,X,Y,PD) 

DGEA0510 

c 

REAL  X,Y(N)  ,PD(1) 

DOBA0520 

c 

DGBA0530 

c 

DOBA0540 

c 

FCNJ  MUST  EVALUATE  PD  IN  BAND  STORAGE  MODE. 

D6BA0550 

c 

THAT  IS,  PD(N*(J-I'fNLC)+I)  IS  THE  PARTIAL  DGBA0560 

c 

DERIVATIVE  OF  YPRIME  (I)  WITH  RESPECT  TO 

DGEA0570 

c 

Y(J) .  NLC  IS  THE  NUMBER  OF  LOWER 

DGBA0580 

c 

CODIAGONALS  FOR  THE  BAND  MATRIX. 

DOBA0590 

c 

FCNJ  MUST  APPEAR  IN  AN  EXTERNAL  STATEMENT  IMDGEA0600 

c 

THE  CALLING  PROGRAM  AND  N,X,Y(1) , . . . ,Y(N) 

DGEA0610 
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X 


H 


Y 


TOL 


HBTH 


MITER 


msr  MOT  BB  ALTBRSD  BY  PCMJ.  D(»MS20 

raur  IS  osxD  omziY  iv  mithi  is  BQmL  to  ocnnosao 

1  OR  -1.  OXBBRMXSX  A  OOMIY  ROOTXn  CAM  1XSA0640 
BB  8QB8TITOTXD.  SBB  RBMMIK  1.  DOAOCSO 

-  IMD8PBMDBMT  VARZBBLB.  (IMPOT  AMD  ODTPOT)  OOBBOSSO 

OM  IWWT.  X  8DPPLZB8  THB  IMITXAL  VAXUB  DOBA0S70 

AMD  IS  USBD  CMOiY  OH  THB  FIRST  CALL.  DCSAOSSO 

ON  ODTPOT,  X  IS  RBKACBD  WITH  THB  CDRRBMT  DOAOSSO 
VALDB  OF  THB  IMDBPBMl»Dffr  VARIABLB  AT  WHICHDGBA0700 

iMracntAnoM  has  bbbm  ccMpLBm).  D(aA07io 

-  IMPOr/OOTPOT.  DSBA0720 

OM  IMPOT,  H  COCTEAIMS  THB  MBXT  STBP  SIZB  IN  DGaA0730 
X.  H  IS  OSBD  ONLY  ON  THB  FIRST  CALL.  DGatA0740 

OH  ODTPOT,  B  CCMTEAIMS  THB  STBP  SIZB  OSBD  D(BA0750 
LAST.  WMBTHBR  SDCCBSSFDLLY  OR  NOT.  DCTA0760 

-  DBPBMDBMT  VARIABLSS,  VBCXOT  OF  LBMCnni  N.  O8BA0770 

(IMPDT  AMD  ODTPOT)  D(»A0780 

OH  IMPDT.  Y(l) .  . . . ,  Y(M)  SUPPLY  INITIAL  DOAOTSO 

VALDBS.  DQBA0800 

OH  ODTPOT,  Ya),...,Y(N)  ARB  RBPLACBD  WITH  DG8A0810 
A  COMPDIBD  VALDB  AT  XBMD.  DnA0820 

-  INPUT  VALDB  OF  X  AT  WHICH  SCODTIOH  IS  DBSIRBD  DGOAOSSO 

MBXT.  IMneB»TIOH  WILL  NORMALLY  00  DOBA0840 

BBYOHD  XBMD  AMD  THB  ROOTIMB  WILL  INTBRPOLATBDGBA0850 
TO  X  ■  XBMD.  DOBAOSeO 

NOTB  THAT  (X*XBMD)  *H  MOST  BB  LESS  THAN  D(SA0870 

ZBRO  (X  AMD  H  AS  SPBCIFIBD  C»  IMPDT)  .  DGOKAOSSO 

•  IMPOT  RBLATIVB  BRROR  BODMD.  TOL  MOST  BB  DCBA0S90 

anSATBR  THAN  ZBRO.  TOL  IS  DSBD  ONLY  ON  THB  DCTAOSOO 
FIRST  CALL  DMLBSS  IMDBX  IS  BQDAL  TO  -1.  DCHUOSIO 
TOL  SaODLD  BB  AT  LBAST  AN  ORDER  OF  D«A0920 

MAONITDDB  LARCOER  TRAM  THB  OMIT  ROUNDOFF  DCBAOSSO 

BUT  OBMXRALLY  HOT  LAROBR  THAN  .001.  DGBA0940 

SIHOLB  STBP  BRROR  BSTIMATBS  DIVTDBD  BY  DCOAOPSO 

YMAXd)  WILL  BB  XBPT  LB8S  THAN  TOL  IN  D(SA09$0 

ROOT-MBAN-SODARB  HORN  (BDCLIDBAN  NORM  DGSAOPTO 

DIVIOBD  BY  SQBrm) )  .  THB  VBCTOR  YMAX  OF  DOBA0980 
WBIGHTS  IS  COHPDTBD  INTBRMALLY  AMD  STORED  DOAOSSO 
IN  WORK  VBCTOR  WK.  IMITIALLY  YMAX  (I)  IS  DOULLOOO 
THB  ABSOLOTB  VALDB  OF  Y(I)  ,  WITH  A  DBRADLT  DOAIOIO 
VALDB  OF  (»B  IF  Y{I)  IS  BQDAL  TO  ZTOO.  D(M»1020 

THBRBAFTBR,  YMAX(Z)  IS  THB  LARGEST  VALUE  DCTAIOBO 
OF  THB  ABSOLUTE  VALUE  OF  Y(I)  SEEN  SO  FAR,  DOBA1040 
OR  THB  INITIAL  VALDB  OF  YMAX  (I)  IF  THAT  IS  DGBAIOSO 
LAROBR.  DCSAIOSO 

•  INPUT  BASIC  METHOD  INDICATOR.  D6BA1070 

USBD  Cm.Y  OR  THB  FIRST  CALL  UNLESS  INDEX  IS  DGBAIOSO 
EQUAL  TO  -1.  DGBAIOSO 

MBTH  -  1,  IMPLIES  THAT  THB  ADAMS  METHOD  IS  D(HCA1100 
TO  BB  USBD.  DCSAlllO 

MBTH  -  2,  IMPLIES  THAT  THB  STIFF  METHODS  OF  DGBA1120 
OBAR,  OR  THB  BACKWARD  DIFFERENTIATION  DOBAliaO 

FORMOLAB  ARB  TO  BB  USBD.  DOBA1140 

•  INPUT  ITERATION  METHOD  INDICATOR.  DOSAllSO 

MITER  -  0.  IMPLIES  THAT  FUNCTIONAL  DOTA1160 

ITERATION  IS  USED.  NO  PARTIAL  DQBA1170 

DBRTVATIVBS  ARE  NEEDED.  A  DUMMY  FCNJ  DGBA1180 

CAN  BE  USBD.  DOBA1190 

MITER  -  1,  IMPLIES  THAT  THB  CHORD  METHOD  DOBA1200 
IS  USED  WITH  AN  ANALYTIC  JACOBIAN.  FOR  DGBA1210 
THIS  METHOD,  THE  USER  SUPPLIES  DOBA1220 

SUBROUTINE  FCNJ.  DGBA1230 
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ISDBX 


imc 

NK 


HXTIR  -  2.  XMPXaXS  tSOa  THE  CHORD  MITBaO  D(»ia240 
IS  U8ID  wrxR  ns  JiiyooBiJua  caxiCDLA!ibd  Dcnunaso 
IHmiDVXiLY  BY  nHITB  DIFFBBBHCBS .  DOAISSO 

A  OOMIY  FCK7  CAR  BB  USBD.  OCaA1270 

MITBR  «  3,  mPLZBS  THAT  THB  CHORD  MBTHCX}  DOAIASO 
IS  DSBD  VITK  THB  JACOBIAN  RBPLACBD  BY  DGSAISSO 
A  DIAflORAL  APPRQXHIATION  BASED  OKI  A  DCBA1300 

DIRBCnOKAXi  DBRIVATIVB.  DCTA1310 

A  DOMMY  FCRJ  CAR  BB  DSBD.  DOBAISSO 

MITBR  •  -1  OR  ~2,  IMPLIBS  DSB  THB  SAME  DGBUa330 

METHOD  AS  ?Ok  MITBR>  1  OR  2,  RBSPBCTIVBLY,DGnUU340 
BUT  USIRO  A  BARDBD  JACOBIAN  MATRIX.  IN  DGRA1350 
THBSB  TNO  CASES  BARDHIDTH  INF(»MATION  DOBA1340 

MOST  BB  PASSED  TO  DOBAR  IHRODCHl  THB  DOBA1370 

comm  BLOCK  D(nA1380 

nniBBIlW  /D8ARD/  RLC.NOC  DOTA1390 

VRBRB  RLC-RDWBR  OF  LONER  CCB>IAQCaiALS  DQBAIAOO 

RDCaRDIBBR  OP  DPPBR  COOIAOONALS  Da8A1410 

INPUT  AND  OOTPUT  PARAMBTBR  DSBD  TO  INDICATE  DaBA1420 
THE  TYPE  OP  CALL  TO  THB  SDBRODTINB.  OR  D(3RA1430 
OUTPUT  IMDBX  IS  RBSBT  TO  0  IP  IRTBORAnOR  D(»A1440 
NAS  SDCCBSSPDL.  OTHERWISE,  THB  VALUE  OF  D6BA1450 
INDBX  IS  tnOlAHQBD.  DQBA146w 

OH  INPUT,  INIHtX  «  1,  IMPLIBS  THAT  THIS  IS  THB  DCaA1470 
FIRST  CALL  FOR  THIS  PROBLEM.  DOBAIASO 

OH  INPUT,  TMDBX  -  0,  IMPLIBS  THAT  THIS  IS  MOT  D<aA1490 
THB  FIRST  CALL  FOR  THIS  mOBLBM.  DOBA1500 

OH  INPUT,  IMDBX  -  -I,  IMPLIBS  THAT  THIS  IS  NOTDOBAISIO 
THB  FIRST  CALL  FOR  THIS  PROBLEM,  AMD  THB  DOBAISSO 
USER  HAS  RBSBT  TOL.  D(aA1530 

ON  INPUT,  INDEX  -  2,  IMPLIBS  THAT  THIS  IS  NOT  DGBA1540 
THB  FIRST  CALL  FOR  THIS  PROBLBM.  IRTBQRATI0HDGBA1550 
IS  TO  COTIMDB  AMD  XBRD  IS  TO  BB  HIT  BXACTLYDOBAISSO 
(NO  INTBRPCMATIOH  IS  DONB)  .  THIS  VALUE  OF  DCTA1570 
IMDBX  ASSOMBS  THAT  XBND  IS  BBYOHD  THB  DQBAISSO 

CURRENT  VALUE  OF  X.  D(«A1590 

OH  IHPUT,  INDBX  -  3,  IMPLIBS  THAT  THIS  IS  MOT  DGBA1600 
THB  FIRST  CALL  FOR  THIS  PROBLBM.  IRTB(»tATIONDGtBA1610 
IS  TO  COMTIMDB  AND  CONTROL  IS  TO  BB  RBTURNEDDGBA1620 
TO  THB  CALLIN8  PROGBUUI  AFTER  OHE  STEP.  XBND  DOTA1630 
IS  lOTORED.  D6BA1640 

IHTBGBR  NORK  VECTOR  OF  LENGTH  N.  USED  ONLY  IF  DGBA1650 
MITBR  ■  1  OR  2  DGBAlSeO 

REAL  NORK  VECTOR  OF  LENGTH  4*N-i-MMBTH-i-NMITER.  DaBAl€70 
THB  VALUE  OF  NMBTH  DEPENDS  ON  THB  VALUE  OF  D(»A1680 
MBTH.  DGBA1690 

IF  MBTH  IS  ^JOAL  TO  1,  OaBA1700 

NMBTH  IS  BQUAL  TO  N*13 .  DGBA1710 

IF  MBTH  IS  EQUAL  TO  2,  DQBA1720 

NMBTH  IS  BQUAL  TO  N*6.  DGBA1730 

THB  VALUB  OF  HMIT8R  DEPENDS  ON  THE  VALUE  OF  DCTA1740 
MITBR.  DaBA1750 

IF  MITBR  IS  B9XAL  TO  1  OR  2,  0(»A1760 

NMITBR  IS  EQPKL  TO  N*(M+1)  DaBA1770 

IF  MITBR  IS  BQUAL  TO  -1  OR  -2,  DOEA1780 

NMITBR  IS  BQUAL  TO  (2*NLC-t-NUC-i-3)  *N  DCTA1790 

NHBRB  NLCoNDimBR  OF  UCmBR  CODIAGONALS  DGBA1800 
NUCaNDMBBR  OF  UPPBR  C(X)IAGONALS  DCHUUSIO 
IF  MITBR  IS  BQUAL  TO  3,  DGBA1830 

NMITBR  IS  BQUAL  TO  N.  D(REA1830 

IF  MITBR  IS  BQUAL  TO  0,  06BA1840 

NMITBR  IS  BQUAL  TO  1.  DGBA1850 
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UK  ID8T  Bimni  DHcrowgiD  bitnibn  socausrvB  oouasso 
au&s  oouas  amunnm.  dqrbisto 

onoR  vMUiMRiR.  {ofjnm)  oaBAisto 

■^ppiTiiry  fpfwpp  pCBIU.890 

ZBR  •  33.  IMWiIBS  TORT  X-^H  VILZ.  BQOAL  X  OR  JCOUISOO 


T8B  HBXr  STBR.  1HZS  GOHDXTZOII  DOBS  NOT 
9QRCB  rat  RCOTZKB  TO  HALT.  BONBVBR,  IT 
DOBS  XNDICRSB  ORB  OF  TWO  CQHDITIOnS. 

TRB  D8BR  MZ<»T  BB  RBQDZRZNO  TOO  MDCH 
ACCORBCT  VIA  TRB  INFOT  PARMOTBR  TOL. 

IB  THIS  CBSB  TRB  OSBR  SHODU)  COBSIDBR 
IBCRBASUIS  rm  VALDB  OF  TOL.  TRB  OTRBR 
GORDinOR  NHICR  MZGHT  GIVB  RISB  TO  THIS 
BRK»  MBSSASB  IS  TRAT  TRB  8TSTBK  OF 
DIFFBRBHTIAL  BQDRTIQBS  BSZNO  SOLVBD 
IS  STIFF  (BZTRBR  IB  SBNBRAL  OR  OVBR 
TRB  SOBZBTBRVAL  OF  TRB  PROBLBM  BBIBG 
solved  at  TRB  TZMB  OF  TRB  BRROR)  .  IN 
TRIS  CASH  TRB  USSR  SBODLD  COBSIDBR 
USINQ  A  BORZBRO  VALDB  FOR  THE  IBPOT 


D(SA1910 

DCTA1930 

DaBA1930 

0(aAI940 

DaBA1950 

D(aA19€0 

DCBUU.970 

DOTA1980 

OGnKA1990 

D(BtA3000 

DSBA2010 

DGSA2020 

OGBA2030 

DaBA2040 

O(»A20S0 

DOBASOSO 

DCTA3070 

DGSA2080 

DQBA3090 

D6BA2100 

D0RA2110 

DORASUO 

DQBA2130 


Step  - 


PRBCISKXr/RARDNARB 


RBQD.  IMSL  RODTINBS 


BARBING  BITR  FIX  BRROR  DCTA3070 

IBR  •  68.  IMFLIBS  TRAT  TRB  BRROR  TBST  DGSA2080 

FAILED.  H  BAS  RBDDCBD  BY  .1  ORB  OR  MORB  DQBA3090 
THOU  AND  TRB  STBF  BAS  TRIBD  AGAIN  D6BA2100 

SDCCBSSFDLLT.  D0RA2110 

IBR  «  67  .  IMFLIBS  TRAT  CORRBCIOR  DOBASUO 

OGHVBRGBBCB  COULD  lOT  BB  ACHIBVBD.  DOBAZISO 

H  BAS  RBDDCBD  BY  .1  OIB  MC»B  TIMBS  AND  DaBA2I40 
TRB  STBF  BAS  TRIBD  AGAIN  SDCCBSSFDLLY .  D(BEA2150 

T»Mmm|T.  ifappe  DCnEA2160 

IBR  •  132.  IMFLIBS  TRB  IBTB6RATION  BAS  D(SA2170 
HALTSD  AFTBR  FAILING  TO  PASS  TRB  BRROR  D(SA2180 
TBST  BVBN  AFTBR  RBDDCING  H  BY  A  FACTOR  DOASISO 
OF  I.OBIO  FROM  US  INITIAL  VALDB.  DOBASSOO 

SBB  RBMARXS.  DOBA2210 

IBR  -  133,  IMFLIBS  TRB  INTBGRATIOR  BAS  DaKA2220 

RALTBD  AFTBR  FAZLIBO  TO  ACRIBVB  DGBA2230 

GCRtRBCTDR  CORVBRGBRCB  BVBN  AFTBR  D(BtA2240 

RBDDCING  R  BY  A  FACTOR  OF  I.OBIO  FROM  DCnA2250 

ITS  INITIAL  VALDB.  SBB  RBMARXS.  DGBA2260 

IBR  »  134,  ZMPLIBS  THAT  AFTBR  SGHB  INITIAL  DOBA2370 
SDCCBSS,  TRB  INTBGRATION  BAS  HALTED  BITHBRDGBA2280 
BY  RBPBATBD  BRROR  TBST  FAILDRBS  OR  BY  DOBA2290 
A  TBST  OH  TOL.  SBB  RBMARXS.  DGBA2300 

IBR  -  135,  IMFLIBS  TRAT  ORB  OF  TRB  INPUT  OOBA2310 
FARAMBTBRS  N,X,H,XBND,TOL,MBTH,MITBR,  OR  D(»A2320 
INDBX  BAS  SFBCIFIBD  INCORRBCTLY.  DGBA2330 

IBR  •  136,  IMFLIBS  THAT  INDBX  HAD  A  VALDB  DCSA2340 
OF  -1  ON  IBFDT,  BDT  TRB  DBSIRBD  CHANGBS  DCaA2350 
OF  FARAMBTBRS  BBRB  ROT  IMPTJaiBNTBD  DGBA2360 

BBCADSB  XBND  BAS  NOT  BEYCBO)  X.  OGBA2370 

IBTBRF(RiAnC»r  TO  X  «  XBND  BAS  PERFORMED.  DGBA2380 
TO  TRY  AGAIN,  SIMPLY  CALL  AGAIN  WITH  DCnA2390 

INDBX  BQDAL  TO  -1  AND  A  NBB  VALDB  FOR  DOBA2400 

XBND.  DGBA2410 

p  -  #  of  integration  steps  taken  * 

DGEA2420 

-  SINGLE  AND  DODBLB/H32  DGBA2430 

-  SINmiB/H36,H48,H60  DCTA2440 

DaBA2450 

DGRCS,DQU3l,DGDtFS,DGRST,LDDATF,LDBZMF,LB0TlB,  DGBA2460 
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NOXATIOH 

RBARICS  1 

2. 

3 


UBXTSTt  OOITZO 

•  ZHTOiaiATIGII  OH  SnCZAL  MOXATIOir  AID 
coMvmnons  is  avazlsbue  m  rm  mkhoal 
maiaoacnoH  c»  Ttafoaaa  hbl  rootthb  ohblp 

1HB  SrriRlIAL  SOBKXmMB  PCK7  IS  OSSD  OEa.Y  1IHDI 
UDDT  PAKAmm  KZTIR  is  igOAL  TO  1  C»  >1.  OTHERinSB, 
A  Down  FOKITIOn  CWI  BB  OSBD.  THB  DUMIY  SOBIiOOTIHB 
fflOOU)  BB  OP  THB  FOLLONZBS  POSH 
SOBRODTnQt  PCHJ  (H.X,Y.PD) 

IWTBQBR  N 

RBAL  Y(II)  .PD(M,N)  ,X 
RBTDRH 


APTBR  THB  INITIAL  CALL,  IP  A  NOmAL  RBTORN  OCCOBKBD 
(IBR-0)  AND  A  NORMAL  CONTIMOATION  IS  DBSIRBD,  SIMPLY 
RBSBT  XBND  AMD  CALL  DQBAR  ABAIN.  ALL  OIHBR 
PARAIOTBRS  WILL  BB  READY  POR  THB  MBXT  CALL.  A  CHAN6B 
OP  PARAMBTBRS  WITH  INDBX  BQOAL  TO  -1  CAN  BB  MADB 
APTBR  BITHBR  A  SOCCBSSPOL  OR  AN  OHSOCCBSSPDL  RBTDRM. 
THB  COMMON  BLOCKS  /DBAND/  AND  /GBAR/  HBBD  TO  BB 
PRBSBRVBD  BBTNBBN  CALLS  TO  D(»AR.  IP  IT  IS  NBCBSSARY 
POR  THB  CdBKSN  BLOCKS  TO  BXIST  IN  THB  CALLING  PROGRAM 


DatA2470 

DOBA2480 

D(aA2490 

D«A2500 

DOBA2510 

DGBA2520 

DGBUa530 

DCTA2S40 

OGBA2550 

DCBA2560 

D(aA2570 

D6BA2580 

D(aA2590 

DCBn2€00 

DGBA3610 

D(SA2£20 

D(SA2€30 

DGaA2640 

1X3A2850 

DG8A2880 

DGBA3670 

DGBA2880 

D(aA2690 

DGatA2700 


THB  POLLOWINS  STATRMBNTS  SHOOU)  BB  INCLDDBD  DCSA2710 

/DBABD/  MLC.NDC  DGBA2720 

/CTAR/  DaMHY(48),SDtBBSY(4)  ,IDCBBfY(38)  008X2730 

WHBRE  DOMDr,  SDCBBIT,  AND  IDDMIY  ARB  VARIABLB  NAMBS  NOT  DOBA2740 
USED  BLSBWHBRB  IN  THB  CALLING  PROCBIAM.  (POR  DOOBLB  DOBA2750 
PRECISION  DOMIY  IS  TYPB  DOOBLB  AMD  SDOMMY  IS  TYPB  RBAL)  DCatA2780 
THB  CHOICB  OP  VALDBS  FOR  MBTH  AMD  MITBR  MAY  RBQOIRB  DGBEA2770 
SOHB  BXPBRIBBMTATKNr,  AMD  ALSO  SCMB  CONSIDBRATTON  OP  DGBA2780 
THB  NATDRB  OP  THB  PROBUtM  AMD  OP  STORAOT  RB(2DIRBHBNTS .  DOBA2790 
THB  PRIMB  CONSIDBRATION  IS  STIPPNBSS.  IP  DOBA2800 

THB  PROBLBM  IS  ROT  STTPP,  THB  BEST  CHOICB  IS  PROBABLY  DOBA2810 
MBTH  «  1  WITH  MITBR  •  0.  IP  THB  PROBLBM  IS  STIPP  TO  A  D(SA2820 
SIGNIFICANT  DB^tBE,  THBN  MBTH  SHODLD  BB  2  AMD  MITBR  DOBA2830 
SHOOU)  BB  1,2, -1, -2  OR  3.  IP  THB  OSBR  HAS  MO  KROWLBDGB  DGBA2840 
OP  THB  IMHBRBNT  TIMB  CONSTANTS  OP  THB  PROBLBM,  WITH  DGSA2850 
WHICH  TO  PREDICT  ITS  STIPPNBSS,  0MB  NAY  TO  DBTBRMINE  DGBA2860 
THIS  IS  TO  TRY  MBTH  >  1  AND  MITBR  >  0  FIRST,  AND  LOOK  D(»A2870 
AT  THB  BEHAVIOR  OF  THB  SOLOTION  CQHPDTBD  AND  THE  STEP  DGBA2880 
SIZES  USBD.  IF  THB  TYPICAL  V«LDBS  OF  H  ARB  MOCH  DCHtA2890 

SMALLER  THAN  THB  SOLDTICNI  BEHAVIOR  WOOLD  SBBM  TO  DCSA2900 

RBQOIRE  (THAT  IS,  MORB  THAN  100  STEPS  ARB  TAKEN  OVER  DOBA2910 
AN  INTERVAL  IN  WHICH  THB  SOLOTICNIS  CHANGE  BY  LESS  DGBA2920 

THAN  ONB  PERCENT)  ,  THEN  THB  PR(SLBM  IS  PROBABLY  STIFF  DGBA2930 
AMD  THE  DEGREE  OF  STIFFNESS  CAN  BB  ESTIMATED  FROM  THE  DQBA2940 
VALDBS  OF  H  USBD  AND  THB  SMOOTHNBSS  OF  THB  SOLOTICNl.  DGBA2950 
IF  THB  DBGRBB  OF  STIFFNESS  IS  ONLY  SLIGHT,  IT  MAY  BE  DOBA2960 
THAT  MBTH»1  IS  MORB  EFFICIENT  THAN  MBTH>2 .  DGBA2970 

EXPBRIMBMTATTON  WOOLD  BB  REQUIRED  TO  DBTBRMINE  THIS.  DGBA2980 
REGARDLESS  OF  MBTH,  THB  LEAST  BPPBCnVE  VALUE  OF  DGBA2990 

MITBR  IS  0,  AMD  THB  MOST  BPPBCTTVB  IS  1,-1, 2, OR  -2.  DGBA3000 
MITBR  •  3  IS  GBNBRALLY  SOHBWHBRB  IN  BBTNBBN.  SINCE  DGBA3010 
THB  STORAGE  REQUIREMENTS  GO  UP  IN  THB  SAME  ORDER  AS  DGBA3020 
BFPBCnVBNBSS,  TRADB-OFF  CONSIDBRATIONS  ARE  DGBA3030 

NBCBSSARY.  FOR  REASONS  OF  ACCURACY  AND  SPBBD,  THB  DOBA3040 

CHOICB  OF  ABS (MITBR) >1  IS  GBNBRALLY  PREFERRED  TO  DGBA30S0 

ABS (MITBR) «2,  UNLESS  THB  SYSTEM  IS  FAIRLY  COMPLICATED  DGBA3060 
(AND  FCNJ  IS  THUS  NOT  PBASIBLB  TO  CODE)  .  THB  DGBA3070 

ACCURACY  OF  THB  FCNJ  CALCULATION  CAM  BE  CHECKED  BY  DGBA3080 
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ooHmsiicM  or  tn  amMoiMi  vzth  tnr  osmuaso  with 
juM(icmnu-3.  ir  m  aMxmvm  mntzx  is  8iaHzricMm.Y 
DliXaOIOftliT  OGKZHMT,  tHBt  IBS  OTTXOI  IIZtSR  -  3  IS 
UXBUt  to  BB  BBBBLT  AS  BffBCTlVB  AS  ABS (MZTBR) -1  OR  2. 
ABD  WILL  SAVB  COBSmBBMWiB  STOAMS  ABD  ROB  TIMB. 

XT  IS  rossiBu.  AID  roxBBrxAix.Y  gom  dbsirablb,  to 

USB  DimtRBBT  VAXJOBS  09  MBXH  ABD  MITBR  IB  DirfBRBirT 
UJBIBTBRVALS  OF  TBB  FRCBLBM.  WOSL  BXMIFIB,  IF  THB 
FROBIBII  IS  BOB-STIFF  IBXTIALI.Y  ABD  STIFF  LATBR, 

MBTH  -  1  ABD  MITBR  •  0  mOBT  BB  SBT  IBITIALLY,  ABD 
MBTH  •  2  ABD  MITBR  -  I  LAIBR. 

5.  THB  INITIAL  VALOB  OF  TBB  STBF  SIBB,  H,  8B00LD  BB 
CHOSBN  COMSIDBRABLY  8MBLLBR  TMMI  THB  AVBRBOB  VALOB 
KIFBCTBD  rC»  THB  PROBUM,  AS  THB  WUtST-COam  MBTHOO 
WITH  WKICK  DOUkR  BBGIBS  IS  MOT  GSBBRALLY  TBB  MOST 
BTFICIBBT  cm.  MOBBVBB.  FOR  THB  FIRST  STBF,  AS  FOR 
BVBRY  STBF.  DGOMR  TB8TS  FOR  THB  POSSIBILITY  THAT 
THB  STBF  SIZB  NAS  TOO  LAROB  TO  PASS  THB  BRRCBl  TBST 
(BASBD  cm  TQL) ,  AMD  IF  SO  ADJOSTS  THB  STBF  SIBB 
DOMM  AOTQMATICALLT.  THIS  DONMNARD  ADJOSTMBHT,  IP 
AMY,  IS  MOIBD  BY  IBR  HAVIMS  TBB  VALOBS  66  OR  67, 

AMD  StBSBQOBMT  ROMS  CM  THB  SAMB  OR  SIMILAR  FROBLBM 
SHCXILD  BB  STARTBO  WITH  AM  AFFROPRIATBLY  SMALLBR 
VALOB  OF  H. 

6.  SOMt  OP  THB  VALOBS  OP  IMTBRBST  LOCATBD  IM  THB 
OCBSKAI  BLOC3C  /(SAR/  ARB 

A.  HOSBD,  THB  STBF  SIBB  H  LAST  OSBD  SOCCSSSFOLLY 

(ODMfY(8)) 

B.  MQOSBD,  TBB  CMtDBR  LAST  OSBD  SOCXXSSFOLLY 

(IIXaBIY(6)) 

C.  MSTBP,  THB  CDMDLATIVB  MOMUBR  OF  STBFS  TAXBM 

(IIXMHfY(7)) 

D.  MFB,  THB  (BBKlLAnVB  MOMOR  OP  FOT  BVALOATIOIIS 

(IDOMMKS) ) 

B.  MOB,  TBB  COMOLATIVB  RBBNER  OF  JACOBIAM 
BVALORTIOHS,  AMD  HBNCB  ALSO  OF  MATRIX  LO 
DBCXMCFOSmasS  dDOMIYO)) 

7.  THB  MC»MAL  OSAOB  OF  DOBAR  MAY  BB  SOMIARIBBD  AS  FOLLOWS 

A.  SBT  THB  INITIAL  VALOBS  IM  Y. 

B.  SBT  M,  X,  K,  TOL,  MBTH,  AK3  MITBR. 

C.  SBT  XBMD  TO  THB  FIRST  OOTFOT  FOUTT,  AMD  IMDBX  TO  1. 

D.  CALL  DGBAR 

B.  BXIT  IF  IBR  IS  CBtBATBR  THAN  128. 

F.  OTRBRNISB,  DO  DBSIRBD  OOTFOl'  OF  Y. 

G.  BXIT  IF  THB  FROBLBM  IS  FIMISBBD. 

H.  OTHBRWISB,  RBSBT  XBMD  TO  THB  MBXT  OOTPOT  POINT,  AMD 

RBTORM  TO  STOP  D. 

8 .  THB  ERROR  WHICH  IS  CGMTROLLBD  BY  WAY  OF  THB  PARAMETER 
TOL  IS  AM  BSTIMATB  OF  THE  LOCAL  TROHCATIOH  BRROR,  THAT 
IS,  THE  ERROR  OOMUTTBU  QM  TAXIMG  A  SIlXnA  STBF  WITH 
THB  MBTHOO,  STARTING  WITH  DATA  RBGARDBD  AS  EXACT.  THIS 
IS  TO  BB  DISTINGUISHED  FROM  THB  GLOBAL  TRUMCATIOM 
FRROR,  WRIC3I  IS  THB  BRROR  IM  AMY  GIVBM  COMPUTED  VALUE 

y/  Y(X)  AS  A  RESULT  OF  THB  LOCAL  TRUMCATIOM  ERRORS 
FROM  ALL  STBFS  TAKBM  TO  OBTAIN  Y(X)  .  THB  LATTER  BRROR 
ACC!DMULATBS  IM  A  MOM-TRIVZAL  NAY  FROM  THB  LOCAL 
ERRORS,  AMD  IS  MBITHBR  BSTIMATBD  NOR  CONTROLLED  BY 
THB  ROOTIMB.  SINCB  IT  IS  USUALLY  THB  (HOBAL  BRROR  THAT 
A  USBR  WANTS  TO  HAVE  UHDBR  CONTROL,  SOB 
BXFBRIMBMTATI(»  MAY  BB  MBCXSSARY  TO  (BT  THB  RIBIT 
VALOB  OF  TOL  TO  ACHIEVE  THB  USBfl3  MBBDS.  IF  THB 
PROBLEM  IS  MATHBHATTCALLY  STABLE,  AMD  THB  METHOD  USED 


D0BA3090 

D0BA3100 

DQBA3110 

OQBA3120 

0(»A3130 

DGHM3140 

DOASISO 

DOBASISO 

D(»A3170 

DGBASISO 

D«A3190 

D(»A3200 

DCHtA3210 

OeaA3220 

DaBA3230 

DGMA3240 

OaBA3250 

D(BA3260 

DOBA3270 

DaBA3280 

DaA3290 

D(»A3300 

D(BA3310 

DaBA3320 

DaBA3330 

D(»A3340 

DGnA3350 

D«A3360 

DGaA3370 

DaA3380 

DQBA3390 

DaBA3400 

DaOBA3410 

DCH(A3420 

DCMA3430 

DaBA3440 

DaBA34S0 

DGMA3460 

DaBA3470 

DaBA3480 

DOBA3490 

OOBA3500 

DaCA3510 

DaA3520 

D(BA3530 

D(BA3S40 

DCTA3550 

DCTA3560 

DGBA3570 

D(BA3580 

D(BA3590 

DGBA3600 

DGBA3610 

D(BA3620 

DGBA3630 

OGBA3640 

DGBA3650 

DOTA3660 

D(BA3670 

DaBA3680 

DGBA3690 

D(»A3700 
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cosmiasrr 

NMtRMmr 


zs  Mnmnxxmx  stmu,  mi  m  olcmm.  man  ga  ji  dcbig7io 
ozvBi  z  SMOLD  vm  SMOOxHur  win  Toi.  i>  A  mmmm  oqbmtzo 
zacHyksmB  *tt— nr  dqbg73o 

zr  tn  Koorzn  Himan  vztr  zzz  vzldbs  Q9  laz,  133,  9CBZ3740 

OR  134.  m  OSIR  8BOOLD  C8ICZ  TO  SB  Zf  TOO  IPCH  OOn37SO 

MsaniKCT  Z8  uzno  notsum.  tb  osn  mkz  «zn  to  Don37co 

SIT  TOZ.  TO  A  XAIOn  VMfli  AM)  COBZIKII.  MOIBIR  DaB3770 

ffOSSZBXA  CaZBB  Of  TWtt  IMOR  COIPZTZOWS  ZS  Al  OOIS37tO 

BWOR  Z«  Tn  COOZMI  09  TH  IITSWM.  fORCTZORS  fCI  lMWk3790 

OR  rcMJ.  zr  MO  moRS  ais  roon.  zr  mux  si  wicisssinr  dohgsoo 
TO  MQHZTOR  ZWTIBMPZIBrS  QOMffZTZlS  BRUtMTSD  ST  THI  DCBA3S10 
ROOrZB.  1BSB  qOiUrrZTZlS  AB  STQRID  ZV  IBI  RORK  VSCTORDCBA3820 

me  AID  imMXSD  BY  sBCzrzc  Muoaom  xa  trs  cciitnif  BLocxDOiAasao 
/QBR/.  zr  zm  zs  132  or  134.  tb  coMromirrs  csoszim  00113140 
TB  error  nsT  rAZums  aui  bi  zonnrzrziD  rROM  xarb  doba3S50 
vRxxiBS  or  TB  QQMrrzrr  dorazsco 

mC(ZOOMMr(ll)-»'Z)/1IK(Z).  rOR  Z-1.....H.  DOIA3870 

OB  CADSI  or  THZS  MAY  B  A  VBY  SBU.  BDT  HOBBRO  DOIA3880 
ZBZnAL  VRZJDB  Of  ABS(T(X))  .  DOBJV3890 

zr  ZBR  ZS  133,  SBVBRAZ.  ROSSZBZLZTZBS  BZZBT.  OBJGSOO 

ZT  BY  B  ZNSTRDCTZVB  TO  TRY  DZrrBRBBT  VAUIBS  Of  1IZTSR.DBA3910 
ALTBRmmVBLY,  TB  OBR  BZOBT  MOSZTOR  SOCCBSSZVB  00183920 

CORRBCTCm  ZTBRATB  OCSRAZRD  ZB  mCCZOOMMYdZ) -fZ) ,  OOBA3930 

Z-1,...,R.  AHOTBR  rOMZBZZiZTY  MZOBT  B  TO  MORZTB  OBA3940 
TB  JAOQBZJUl  mXRZZ,  Zf  OB  ZS  DSBO,  STORBO,  BY  OOBA3950 

caLom,  ZB  mc(zoiaaar(iouz),  roR  z-i,....Bni  zr  oobbsso 

AB(MZTBR)  ZS  BOORL  TO  1  OR  2,  OR  ftBl  Z.l . B  Zf  OBA3970 

MZTBR  ZS  BQIIAL  TO  3.  OORA3980 

OOBA3990 

•  1984  BY  ZMSL,  ZBC.  AU.  RZOBTS  RBSBRVBO.  OOBA4000 

OOBA4010 

•  ZMSZ.  NARRABTS  ORLY  THAT  ZMSL  TBSTZBO  HAS  BBBB  OBA4030 

AFPLZBO  TO  THZS  OC»B.  HD  OTHBR  BARRABTT,  OOBM030 
BXPBSSBO  cm  ZMRLZBO,  ZS  APFLZCABLB.  00184040 


OC«A4070 

SOBROOTZB  OOBAR  (H,rCH,rCHJ,Z,H,Y,XBBO,TOL,MBTH,MXTBR,ZBOBX,  OOBA4080 

1  ZBK,1fK,ZBR,>tm?)  + 

SPBCZrZC»TZaRS  FB  AROOBBTS  OBA4100 

ZBTBBR  B,MBTH,MZTBR,ZBDBZ,Z1IK(1)  ,IBR,Bt«p 

OOOBLB  PRSeZSKW  X,H, Y(B)  ,ZBllD,TOL,mC(l)  OC8A4120 

SPBCZrZCATZGHS  FOR  LOCAL  VARIABLES  OOBA4130 
IBTBQSR  HBRROR,HSAVBl,BSAVB2.RPH,RY,BC,MrC,KrLAS,  DOBA4140 

1  J8TART,BSQ,HgDSBO,RSTBr,HFB,HJB,I,BO,BHCOT,KGO,OOBA4150 

2  JBR,KBR,HB,HBQOIL,IIRaBfY(21)  ,HLC,MDC  OBA41SO 

RBAL  SOaiRIY(4)  OGBA4170 

OOOBLB  PRBCISKm  T,HH,HKIB,HMAX,BPSC,ORODHO,BPSJ.HOSBO,TOOTP,  00884180 

1  AYZ,0,0H,SBPS,DaMMYf39)  OBA4190 

BXXBRMAL  PCB.FCBJ  OGBA4200 

rTBOT  /BAHO/  HLC.BDC  OGBA4210 

GOMMOH  /OBAR/  T,HH,HMZB,HMAX,BPSC,URODHO,BPSJ,HOSBO,OaHMY,  OGBA4220 

1  TOaTP,SODIBlY,IK;,MFC,KFLAG,JSTART,HSQ,HQDSBO,  OOBA4230 

2  HSTBP,BFB,HJB,HPH,HBRR0R,NSAVB1,HSAVB2,NBC20IL,  DOBA4240 

3  BY.ZOOMCY.HO.BBCOT  OGBA4250 

DATA  SBPS/Z3410000000000000/  DBA4260 

FIRST  EZBCDTABLE  STATBMBBT  DOBA4270 

ZF  (MZTBR. GB.O)  BLC  -  -1  DGBA4280 

XBR  -  0  DOBA4290 

OB  .  0  DOBA4300 

ORCOHD  •  SBPS  DBA4310 

CGRPOTB  WORK  VBCTIOR  INDICIBS  DGBA4320 


DATA 
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1 

I' 

■natCK  -  ■ 

D(aA4330 

1 

MKvn  -  iiino»fii 

D(aA4340 

mtetm  -  iisiraii<i« 

DGntA4350 

1 

m  m  IBAVI2+W 

OCSA43CO 

1 

XV  (MtH.IQ.l)  nOOXL  •  liX4>13*ll 

D(aA4370 

t 

IV  (iBTH.iQ.a)  aoKioxx.  - 

DCaA4380 

mw  m  nQoxL  *  H 

DOBA4390 

IV  (1IXTIR.BQ.0.(».IIITBR.BQ.3)  MW  -  MBQOIL 

DOBA4400 

MVC  -  lO^MRM-t-IABS  (MmR) 

OGBA4410 

c 

CHSCaC  TOR  IMCCMOtSCT 

INPUT  PARAMETERS 

D(aA4420 

c 

DWA4430 

IV  (mriR.LT.-a.OR.icxTiR.aT.i)  oo  to  as 

OOBAA440 

IV  (linH.MB.l.AHD.MrrH.MB.2)  00  TO  85 

DOBA4450 

IV  (TOL.LB.O.DO)  00  TO  85 

OOBA44SO 

XV  (H.LB.O)  00  TO  85 

DaBA4470 

IV  ((X-XSVD)*H.OI.O.DO)  00  TO  85 

D(aA4480 

IV  (IMDBX.VQ.O)  00  TO  10 

DOBA4490 

IV  (IMDMX.BQ.S)  00  TO  15 

D(»A4500 

IV  (IMDBX.XQ.-l)  00  TO  20 

OOBA4510 

IF  (I1ID8X.BQ.3)  00  TO  25 

DWA4S20 

IF  (IMDXZ.MB.l)  00  TO  85 

IX3BA4530 

c 

IP  IMITIAIi  VAUnS  OF  WAX  OTHHR  THAN 

DOTA4540 

c 

THOSB  SBT  BBLOW  ARB  DBSIRBD, 

THEY 

OOBA4550 

c 

SHCOLD  BB  SBT  HBRB.  ALI.  WAX(I) 

DOBA4560 

c 

Mrar  BB  posmvB. 

IF  VALDES 

FOR 

DGBA4570 

c 

HMXM  OR  WAX,  THB 

BODHDS  on 

DOBA4580 

c 

DABS(HR) ,  OTHBR  THAM  TROSB  BBXOW 

DOBA4590 

c 

BRB  OBSIRBD,  TOBY 

SHOULD  BB 

SBT 

OOBA4600 

c 

BBXOW. 

DOBA4610 

00  5  I«1,N 

D(3BAA620 

nCd)  >OABS(Y(I)) 

DGBA4630 

XV  (WKd)  .BQ.o.oo)  mcd)  .  i.oo 

D(SA4640 

KKdnr+i)  •  Yd) 

DGXA4650 

5  OOHTIMOB 

DOBA4660 

HC  »  R 

DOBA4870 

M 

a 

DGBA4680 

HH  «  H 

DGBA4690 

IV  ((T+HH)  .BQ.T)  RSR  >  33 

DOBA4700 

HMIM  >  DABS(H) 

DOBA4710 

miAX  •  DABS (X-XBMD) *10.00 

DOBA4720 

BPSC  «  TOL 

DOBA4730 

JSTART  -  0 

DGBA4740 

NO  >  N 

VK3MM150 

HSQ  >  HO*NO 

DOTA4760 

BVSJ  -  DS(^(ORODND) 

DOBA4770 

NRCDT  -  0 

DOBA4780 

oaiafY(2)  «  1.000 

OOBA4790 

DOMMYda)  -  1.000 

DOBA4800 

00  TO  30 

DGtBA4810 

c 

TOOTP  IS  THB  PRBVIOOS  VALUE  OF  XBND 

DOBA4820 

c 

FOR  USB  IN  WAX. 

DOBA4e30 

10  HMAX  «  DABS (XBND-TOUTP) *10.00 

DOBA4840 

OO  TO  45 

DGBA4850 

c 

DOBA4860 

15  HMAX  >  DABS  (XBND-TOUTP)  *10.00 

DWA4870 

IF  ((T-XBND)*HH.OE.O.DO)  GO  TO  95 

DOBA4880 

GO  TO  50 

DCSA4890 

c 

DOBA4900 

20  IF  ( (T-XBMD)*HH.OB.O.DO)  00  TO  90 

DOBA4910 

JSTART  -  -1 

D(SA4920 

NC  >  N 

DOBA4930 

BPSC  «  TOL 

DOBA4940 

56 

f"  ■ 

c 

OGEA4950 

35 

IF  ( (T+HH)  .IQ.T)  KIR  .  33 

DaEA4960 

writ* (*,*), '•rror  cod*  ■  ',k*r 

♦ 

c 

DOEA4970 

30 

IW  >  HO 

D(»A4980 

•t^  ■  stop  +  1 

+ 

writ* (*,*) 'stop  >  ',at*p 

+ 

call  natST  (FCH.FC3KJ,inC(Hy-i>l)  .HK.HKCRBSROR+I)  .mCCNSAVSl+l)  , 

DfflUl4990 

1  NKdisjivBa-t-i)  ,mc(NPN-fi)  .mcdngoiL+i)  ,iinc,NN,8tttp) 

+ 

c 

DQEA5010 

RQO  •  l-KFLAO 

D(NUV5020 

GOTO  (35,55,70,80),  KOO 

OGOASOSO 

c 

ICFLAG  .  0,  -1,  -2,  -3 

DOEA5040 

35 

coRnnoB 

D6BA5050 

c 

nORMAL  RBTDRN  FROM  IRTBORATOR.  THE 

DG«A50€0 

c 

NIKSiTS  YWam  ARB  DPDATBD.  IF 

OCaA5070 

c 

DIFFBRBirr  VALDBS  ARB  DBSIRBD,  THEY 

D(»A5080 

c 

SmOLD  BB  SET  HERB.  A  TEST  IS  MADE 

DOBA5090 

c 

VCSL  TOL  BBIHG  TOO  SMALL  FOR  THE 

D6BA5100 

c 

MACHIHB  PRECISION.  ANY  OTHER  TESTS 

DGBASllO 

c 

OR  CALCULATIONS  THAT  ARB  REQUIRED 

DGEAS120 

c 

AFTER  EVERY  STEP  SMOOLD  BB 

DQEA5130 

c 

INSERTED  HERB.  IP  INDEX  «  3,  Y  IS 

D6EAS140 

c 

SET  TO  THE  CDRRENT  SOLUTION  ON 

DGBA5150 

c 

RETURN.  IF  INDEX  >  2,  HH  IS 

D(3EA5160 

c 

CONTROLLED  TO  HIT  XEND  (WITHIN 

DCTA5170 

c 

ROONDOFF  ZKROXL)  ,  AND  THEN  THE 

0(»A5180 

c 

CDRRENT  SOLDTION  IS  POT  IN  Y  (Xt 

D«A5190 

c 

RETURN.  F(»l  ANY  OTHER  VALOB  OF 

D(SA5200 

c 

INDEX,  CONTROL  RETORHS  TO  THE 

D6EA5210 

c 

INTSGRAKXl  UNLESS  XEND  HAS  BEEN 

DOEAS220 

c 

REACHED.  THEN  INTERPOLATED  VALDES 

DCTA5230 

c 

OF  THE  SOLUTION  ARB  COMPUTED  AND 

DQEA5240 

c 

STORED  IN  Y  ON  RBTDRN. 

DGSA5250 

c 

IF  INTERPOLATION  IS  NOT 

DOBA5260 

c 

DBSIRBD,  THE  CALL  TO  DGRIN  SHOULD 

DGBA5270 

c 

BE  REMOVED  AND  CONTROL  TRANSFERRED 

DGBA5280 

c 

TO  STATEMENT  95  INSTEAD  OF  105. 

DOEA5290 

D  «  O.DO 

DGEA5300 

DO  40  I«1,N 

D(SA5310 

AYI  -  DABS(NK(NY-l-I)  ) 

DGBA5320 

mc(i)  -  ixiAXKincd)  ,AYi) 

D6EA5330 

40 

D  «  D+(AYI/NK(I) )**2 

DQBA5340 

D  -  D* (URODND/TOL) **2 

D6EA5350 

DN  «  N 

DGEAS360 

IF  (D.GT.DN)  GO  TO  75 

OGBA5370 

IF  (INDEX. EQ. 3)  GO  TO  95 

DGBA5380 

IF  (INDEX. EQ. 2)  GO  TO  50 

D6EA5390 

45 

IF  ( (T-XEND)*HH.LT.0.D0)  GO  TO  25 

DGBA5400 

NN  >  NO 

DOBA5410 

CALL  DGRIN  (XBND,HK(NY-l-l)  ,NN,Y) 

DGBA5420 

X  -  XEND 

D6BA5430 

GO  TO  105 

DGBA5440 

50 

IF  (( (T-t-HH) -XEND)  *HH. LB. 0.00)  GO  TO  25 

DGBA5450 

IF  (DABS(T-XBHD)  .LB.UROUND*DMAX1(10.DO*DABS(T)  ,HMAX) )  GO  TO  95 

DGBAS460 

IF  ((T-XEND)*HH.6B.0.D0)  GO  TO  95 

DGBA5470 

HH  >  (XBND-T)*(1.D0-4.D0*URCX]ND) 

DGBA5480 

JSTART  .  -1 

DGBAS490 

GO  TO  25 

DOBA5500 

c 

ON  AN  ERROR  RETURN  FRCBl  INTEGRATOR, 

DGBAS510 

c 

AN  UOIBDIATE  RETURN  OCCURS  IP 

DGBA5520 

c 

KFLAO  -  -2,  AND  RECOVERY  ATTEMPTS 

DGBAS530 
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ASB  MMDB  OXRnanSB.  TO  RBCX3VBR,  HH 
MB)  am  ARB  RBDOCBD  BY  A  FACTOR 
OP  .1  UP  TO  10  TUBS  BBF(aB  01VIH6 
UP. 


55  JBR  -  €5 

60  IP  (NHCUT.BQ.IO)  GO  TO  65 
NHCUT  -  RHCUr-fl 
am  >  aiiN*.ioo 
ai  «  HH*.1D0 
JSTART  ■  -1 
GO  TO  25 

65  IP  (JBR.BQ.66)  JBR  -  132 
IP  (JBR.BQ.67)  JBR  >  133 


DaBA5540 

D(3BA5550 

D(aA5560 

D(aA5570 

0(»A5580 

DaBA5590 

DCaA5600 

DOBA5610 

D(»A5620 

DOBA5630 

OGBA5640 

DCTA5650 

DGBA5660 

D(aA5670 


C 

c 

c 

c 

c 


GO  TO 

95 

DGBA5680 

DGRA5690 

70 

JBR  - 

134 

DGBA5700 

GO  TO 

95 

DGBA5710 

DGBA5720 

75 

JBR  - 

134 

DOBA5730 

KPLAG 

-  -2 

DaBA5740 

GO  TO 

95 

DGBA5750 

DGBA5760 

80 

JBR  - 

67 

DGBA5770 

GO  TO 

60 

D6BA5780 

DOBA5790 

85 

JBR  « 

135 

DGBA5800 

GO  TO 

110 

OOBA5810 

DaBA5820 

90 

JBR  • 

136 

DCaA5830 

RH  >  NO 

CALL  DGRIN  (XEND^NKiNY-fl)  ,NN,Y) 
X  >  XBRD 
GO  TO  110 


95 


INDEX  «  KPLAG 


X  -  T 

DO  100  I«1,N 
Y(I)  -  WK(NY-i-I) 

IF  (JBR.LT.128) 

TOUTP  «  X 
IF  (KPLAG. BQ.O) 

IF  (KPLAG. NE.O) 

IBR  -  MAX0(KBR,JBR) 

CONTIHUB 

IP  (KBR.NB.O. AND. JBR.LT.128)  CALL  UBRTST  (KBR,6HDGBAR  ) 
IP  (JBR.NB.O)  CALL  UBRTST  (JBR,6HDGBAR  ) 

9005  RETURN 
BND 


100 

105 


110 

9000 


H 

H 


HUSBD 

HH 


OGBA5840 

DOBA5850 

DGBA5860 

DaRA5870 

DOTA5880 

DGBA5890 

DOBA5900 

D6BA5910 

DGBA5920 

DOBA5930 

DGBA5940 

DGBA5950 

DGEA5960 

DGBA5970 

DGBA5980 

DGEA5990 

DGEA6000 

OGEA6010 
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ZMSL  ROOTINB  KAMB  •  DGRST  DCHISOOIO 

DGRS0020 

-aodified  to  print  sheath  and  diagnostic  output  to  files  "sheatha.dat" 


and  "diag.dat” 
COMPOTBR 

LA.TBST  RBVISIOM 

PORPOSB 

PRBCISIQN/HARDNARB 


-  ZBM/DOOBLE  DGRS0050 

DGRSOOeO 

-  JUNB  1,  1982  DaRS0070 

DORS0080 

-  NDCLBUS  CALLBD  ONLY  BY  IMSL  SUBROOTZNE  DGBAR  DGRS0090 

DORSOlOO 

-  SINGLB  AND  DOOBLB/H32  DGRSOllO 

-  SINQX.B/H36,H48,H60  DaRS0120 


-  JUNB  1,  1982 


RE(P.  ZMSL  ROUTZNBS  -  OGRCS , DORPS , LUDATF , LDBLMF , LEQTIB , UBRTST , 

DGBTZO 


NOTATZON 


COPYRZGHT 


WARRANTY 


DGRSOllO 

DORS0120 

DGRS0130 

DGRS0140 

DORSOISO 

DGRS0160 

DGRS0170 

DGRS0180 


SUBROUTINE  DGRST 


ZNTEGBR 

DOUBLE  PRSCZSZON 


INTEGER 


REAL 

DOUBLE  PRECISION 


EXTERNAL 
COMMON  /DBAND/ 
COMMON  /GEAR/ 


Open (units8 , file: 
open  <units9 , file: 
KFLA6  «  0 
TOrj)  -  T 


-  INFORMATION  ON  SPECIAL  NOTATZON  AND  DGRS0170 

CONVENTZONS  ZS  AVAZLABLE  ZN  THE  MANUAL  DGRS0180 

ZNTRODUCTZON  OR  THROUGH  ZMSL  ROUTZNE  UHELP  DGRS0190 

DORS0200 

-  1982  BY  ZMSL,  ZNC.  ALL  RZGHTS  RESERVED.  DGRS0210 

DGRS0220 

-  ZMSL  WARRANTS  ONLY  THAT  ZMSL  TESTZNG  HAS  BEEN  DGRS0230 

APPLZED  TO  THIS  CODE.  NO  OTHER  WARRANTY,  DGRS0240 

EXPRESSED  OR  ZMPLZED,  ZS  APPLZCABLE.  DORS0250 

DORS02S0 

. DGRS0270 

DGRS0280 

(FCN,FCNJ,Y,YMAX, ERROR, SAVE1,SAVE2,PW,EQUZL,  DGRS0290 

ZPZV,NO,Step)  * 

SPECZFZCATZONS  FOR  ARGUMENTS  DGRS0310 

ZPZV(1),N0  DGRS0320 

Y(N0,1) ,YMAX(1) ,SRROR(l) , SAVEl ( 1 ) , SAVB2 (1) ,  DGRS0330 

PW(l)  ,EQUIL(l)  ,epriiiie,epriine(2) 

SPBCZFICATZONS  FOR  LOCAL  VARZABLES  DGRS0350 

N,MF,KFLAG,JSTART, HOUSED, NSTEP,NFE,NJE,NSQ,  DGRS0360 

Z,MBTH,MZTBR,nQ,L,ZDOUB,MFOLD,NOLD,ZRET,HBO,  DGRS0370 

MZO,ZHBVAL.MAXDBR,ZJiAX,ZRBDO,  J,NSTEPJ,  Jl,  J2,  DGRS0380 

M,  ZER, NBWQ, NPW,NBRROR, NSAVBl , NSAVB2 , NBQUZL, NY,  DGRS0390 

MZTBR1,ZDUMMY(2) ,NLC,NUC,MWK, JER  DGRS0400 

TQ(4)  D6RS0410 

T,H,HMZN,HMAX,EPS,UROUND,HUSED,EL(13) ,OLDLO,  DGRS0420 

TOLD, RMAX,RC, CRATE, EPSOLD, HOLD, FN,EDN,E,EUP,  DORS0430 

BND, RH, R1 , CON, R,HLO , RO, D, PHLO , PR3 , D1 , BNQ3 , BNQ2 , DGRS0440 
PR2,PR1,ENQ1, EPS J, DUMMY, tcum  + 

FCN,FCNJ  DGRS0460 

NLC,NUC  DGRS0470 

T,H,HMZN,HMAX,EPS,UROUND,EPSJ,HUSED,  DGRS0480 

BL, OLDLO , TOLD , RMAX, RC, CRATE , EPSOLD , HOLD , FN,  DORS0490 

EDN,E,EUP, BND, RH,R1,R,HL0,R0,D, PHLO, PR3,D1,  DGRS0500 

BNQ3,ENQ2,PR2,PR1,ENQ1, DUMMY, TQ,  DGRS0510 

N,  MF , KFLAG , JSTART, NSQ, NQUSED , NSTEP , NFE , NJE ,  DGRS052  0 

NPW,NERROR,NSAVEl,NSAVE2, NBQUZL, NY,  DGRS0530 

Z,METH,MZTER,NQ,L,ZDOUB,MFOLD,NOLD,ZRET,HEO,  DGRS0540 

MZO,ZWEVAL,MAXDER,LMAX,ZREDO, J,NSTEPJ, Jl, J2,  DGRS0550 

M,NEWQ,ZDUMMY  DGRS0560 

FIRST  EXECUTABLE  STATEMENT  DGRS0570 

•'sheatha.dat' ,statuss' unknown' ) 

: '  diag .  dat ' ,  statuss '  \inknown' ) 


THZS  ROUTZNE  PERFORMS  ONE  STEP  OF 


DGRS0580 

DGRS0590 

DGRS0600 
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5 


10 


IP 

IP 


(JSTART.GT.O) 

(JSTART.MB.O) 


GO  TO  50 
GO  TO  10 


THB  niTBGRATIOa  OP  AN  INITIAL  D(»tS0610 

VALOB  PR£»LBM  POR  A  STSTBM  OP  D6RS0620 

ORDINARY  DIPPBRKWriAL  BQOATICNIS.  D«IS0630 

DaRS0640 

ocatsosso 

ON  THB  FIRST  CALL,  THE  ORDBR  IS  SET  DGRS0660 
TO  1  AND  THE  INITIAL  YDOT  IS  DGRS0670 

CALCULATED.  RMAX  IS  THB  MAXIMUM  DGRS0680 
RATIO  BY  WHICH  H  CAN  BE  INOtBASED  DGRS0690 
IN  A  SINGLE  STEP.  IT  IS  INITIALLY  DGRS0700 
1.B4  TO  COMPENSATE  FOR  THB  SMALL  ZX»S0710 
INITIAL  H,  BUT  THEN  IS  NORMALLY  D6RS0720 
BOUAL  TO  10.  IP  A  FAILURE  OCCURS  DORS0730 
(IN  CX>RRBCTOR  CONVERGENCE  OR  ERROR  DGRS0740 
TBST) ,  RMAX  IS  SET  AT  2  POR  THB  DGRS0750 
NEXT  INCREASE.  OGRS0760 


CALL  FCN  (N.T.Y.SAVBl 

DO  5  I«1,N 

Y(I,2)  -  H*SAVB1(I) 

MBTH  «  MP/10 

MITER  >  MF-10*MBTH 

NQ  «  1 

L  >  2 

IDOUB  >  3 

RMAX  «  1.D4 

RC  «  O.DO 

CRATE  a  l.DO 

HOLD  -  H 

MPOLD  •  HP 

NSTBP  -  0 

NSTEPJ  a  0 

NFB  a  1 

NJE  a  0 

IRBT  a  3 

GO  TO  15 


,  eprime ,  epriiBe2 ) 


+ 


DGRS0780 

DGRS0790 

DGRS0800 

OGtRS0810 

DGRS0820 

DORS0830 

DGRS0840 

DGRS0850 

DGRS0860 

DGRS0870 

D6RS0880 

D(H(S0890 

D6RS0900 

D(»S0910 

DaiS0920 

OQRS0930 

DQRS0940 

DGRS0950 

IP  THE  CALLER  HAS  CHANGED  METH,  DCtRS0960 

DGRCS  IS  CALLED  TO  SET  THE  DGRS0970 

COEFFICIENTS  OP  THB  METHOD.  IF  THE  DGRS0980 
CALLER  HAS  CHANGED  N,  EPS,  OR  DGRS0990 

MBTH,  THB  CCEfSTANTS  E,  EDN,  BUP,  D6RS1000 
AND  BND  MUST  BE  RESET.  B  IS  A  DGRSlOlO 

CCEIPARISON  POR  ERRORS  OF  THB  DGRS1020 

CURRENT  ORDER  NQ.  EUP  IS  TO  TBST  DGRS1030 
POR  INCREASING  THB  ORDER,  EDN  FOR  DGRS1040 
DECREASING  THE  ORDBR.  BND  IS  USED  D6RS1050 
TO  TBST  FOR  C0RVBRGBNC:E  OF  THE  DGRS1060 
CORRECTOR  ITERATES.  IF  THE  CALLER  DGRS1070 
HAS  CHANGED  H,  Y  MUST  BE  RESCALED.  DGRS1080 
IP  H  OR  MBTH  HAS  BEEN  CHANGED,  DGRS1090 

IDOUB  IS  RESET  TO  L  1  TO  PREVENT  DGRSllOO 
FURTHER  CHANGES  IN  H  FOR  THAT  MANY  DGRSlllO 
STEPS.  DGRS1120 


IF  (MF.EQ.MFOLD)  GO  TO  25 

MEO  a  MBTH 

MIO  a  MITER 

MBTH  a  MF/10 

MITBR  a  MF-10*MBTH 

MPOLD  a  MF 

IF  (MITER. NB. MIO)  IWEVAL  a  MITER 
IF  (MBTH. EQ. MEO)  GO  TO  25 
IDOUB  a  L+l 
IRET  a  1 


DGRS1130 

DGRS1140 

DGRS1150 

DGRS1160 

D6RS1170 

DGRS1180 

DGRS1190 

DGRS1200 

DGRS1210 

DGRS1220 
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15  CALI.  DOKCS  (lOmi.NQ,  BL.TQ.UAXOIR) 

UnX  -  MAXOBR-i-l 
RC  >  RC*1L(1)/01.DLO 
OLDLO  -  BL(1) 

20  FN  •  N 

BDH  ■  FM*(TQ(1)*BPS)**2 
B  -  ra*(TQ(2)*BPS)**2 
BDP  -  PN*(TQ{3)*BPS)**2 
BUD  .  PN*(TQ(4)*BPS)**2 
BPSOLD  -  BPS 
NCnj)  «  N 

GO  TO  (30,35.50) .  ZRBT 

25  IF  ( (BPS. BQ. BPSOLD)  .AND. (N.BQ. HOLD) )  GO  TO  30 
IF  (N.EQ.NOLD)  IWBVAL  «  MITBR 
IRBT  a  1 
GO  TO  20 

30  IF  (H.BQ.HOLD)  GO  TO  50 
RH  -  H/HOLD 
H  -  HOLD 
IRBDO  >  3 
GO  TO  40 

35  RH  -  DMAX1(RH,HMIN/DABS(H) ) 

40  RH  >  DMIN1(RH,HMAX/DABS(H)  ,RMAX) 

R1  a  l.DO 
DO  45  J-2,L 
R1  >  R1*RH 
DO  45  Ial,N 

45  Y(I,J)  «  Y(I,J)*R1 
H  -  H*RH 
RC  a  RC*RH 
IDOOB  a  L+l 

IF  (IRBDO. BQ.O)  GO  TO  285 

THIS  SBCTION  CGMPDTBS  THB  PRBDICTTBD 
VALDES  BY  BFFBCTIVBLY  MULTIPLYING 
THB  Y  ARRAY  BY  THB  PASCAL  TRIAMGLB 
MATRIX.  RC  IS  THB  RATIO  OF  NBH  TO 
OLD  VALDBS  OF  THB  COEFFICIBHT 
H*BL(1) .  NHBN  RC  DIFFBRS  FROM  1  BY 
IX}RB  TRAN  30  PBRCBNT,  OR  THB 
CAI.T.BR  HAS  CHANGED  MITBR,  INEVAL 
IS  SET  TO  MITBR  TO  FORCE  THE 
PARTIALS  TO  BE  UPDATED,  IF 
PARTIALS  ARE  USED.  IN  iUnr  CASE, 

THB  PARTIALS  ARE  UPDATED  AT  LEAST 
EVERY  20-TH  STEP. 

50  IF  (DABS(RC-l.DO) .6T.0.3D0)  IWBVAL  a  MITBR 
IF  (NSTEP.GE.NSTEPJ+20)  IWBVAL  a  MITBR 
T  a  T+H 
DO  55  Jlal.NQ 
DO  55  J2aJl,NQ 
J  a  (NQ+Jl) -J2 
DO  55  lal.N 

55  Y(I,J)  a  Y(I, J)+Y(I, J+1) 

UP  TO  3  CORRECTOR  ITERATIONS  ARB 
TAKEN.  A  CONVERGENCE  TEST  IS  MADE 
ON  THB  R.M.S.  NORM  OF  EACH 
CORRECTION,  USING  BND,  WHICH  IS 
DEPENDENT  ON  BPS.  THB  SUM  OF  THE 
CORRECTIONS  IS  ACCUMULATED  IN  THE 
VECTOR  ERROR(I)  .  THB  Y  ARRAY  IS 
NOT  ALTERED  IN  THE  CORRECTOR  LOOP. 
THB  UPDATED  Y  VECTOR  IS  STORED 


D«tS1230 

D(atS1240 

DGBiS1250 

DGBIS1280 

O(MtS1270 

DCBtS1280 

DCntS1290 

DORS1300 

DGRS1310 

DORS1330 

DQRS1330 

DGRS1340 

DGRSISSO 

DGRS1360 

DratS1370 

DGRS1380 

DGRS1390 

DGRS1400 

DORS1410 

D(3RS1420 

D(MtS1430 

DCHIS1440 

D(HIS1450 

D6RS1460 

OWXS1470 

D(«S1480 

DGRS1490 

D(ntS1500 

OORSISIO 

lXniS1520 

D(HIS1530 

DGRS1540 

DGRS15S0 

DORS1S60 

DORS1570 

DGRS1S80 

D(HtSlS90 

D6RS1600 

DGKS1610 

DGRS1620 

DORS1630 

DGRS1640 

DGRS1650 

DGRS1660 

D(HtS1870 

D6RS1680 

DGRS1690 

DGRSX700 

DGRS1710 

DGRS1720 

DGRS1730 

DaRS1740 

D^1750 

DGRS1760 

DGRS1770 

DGRS1780 

DGRS1790 

DGRS1800 

DGRS1810 

D6RS1820 

DGRS1830 

DGRS1840 
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C 


TIMPaUlRILY  IH  8AVB1. 

CO  DO  65  I-l.H 
65  UnORd)  -  0.00 

M  •  0 

CALL  rCN  (N,T,Y,SAVB2,epriflM.0priaM2) 

NR  «  NR-i-1 

IF  (IWSVAL.LB.O)  GO  TO  95 

IF  INDICATED,  THB  MATRIX  P  -  I  • 
H*BL(1)*J  IS  RBBVALOATBD  BSFORB 
STARTING  THB  CORRBCTOR  ITBRATION. 
mVAL  IS  SBT  TO  0  AS  AN  INDICA‘lX>R 
THAT  THIS  HAS  BBBN  DONB.  IF  HTTBR 
»  1  OR  2.  P  IS  COMPOTBD  AMD 
PROCBSSBD  IN  PSET.  IF  MITBR  -  3, 
THB  MATRIX  USBD  IS  P  •  I  • 
R«BL(1)*D,  NHBRB  D  IS  A  DIAGONAL 
MATRIX. 

IWBVAL  >  0 
RC  «  l.DO 
NJB  -  NJB+l 
NSTBPJ  -  NSTBP 
GOTO  (75,70,80),  MITBR 
70  NR  >  NR+N 
75  CCRI  a  -H*BL(1) 

MITBRl  «  MITBR 

CALL  D6RPS  (FCH,  FCNJ,  Y, NO ,  CON, MITBRl ,  YMAX,  SAVBl ,  SAVB2 ,  PW,  BQOIL, 

1  IPIV,IBR) 

IF  (IBR.NB.O)  GO  TO  155 
GO  TO  125 
80  R  «  BL(1)*.1D0 
DO  85  I«1,N 

85  PW(I)  «  Y(I,1)+R*(H*SAVB2(1)-Y(I,2)) 

CALL  FCN  (N,T,PN, SAVBl, •prine,epriBM2) 


DGRSliBO 

OC»81860 

DORS1870 

DCKS1680 

♦ 

0(3851900 

0(atS1910 

0^1920 

DORS1930 

DGRS1940 

0(atS1950 

0^1960 

DGHtS1970 

OORS1980 

D(3RS1990 

DaRS2000 

DaRS2010 

DaRS2020 

DGHtS2030 

D(»S2040 

DGBtS2050 

0(»tS2060 

0^2070 

DCMtS2080 

DGRS2090 

D(nS2100 

D(»S2110 

D(»S2120 

DGRS2130 

00182140 

DaRS2l50 

DGRS2160 


NR  •  NR+l  DORS2180 

HLO  «  H*BL(1)  DOtS2190 

DO  90  I«1,R  DORS2200 

RO  «  H*SAVB2(I) -Y(I,2)  D(3R52210 

PW(I)  -  1.00  DGRS2220 

D  -  .1D0*R0-H*(SAVB1(I) -SAVB2(I))  D(3RS2230 

SAVBl (I)  =  O.DO  DORS2240 

IF  (DABS(RO)  .LT.UROORD*YMAX{I} )  60  TO  90  OGRS2250 

IF  (DABS  (D)  .EQ. O.DO)  GO  TO  155  D6RS2260 

PN(I)  >  .100*RO/D  DGRS2270 

SAVBl (I)  «  PH(I)*RO  DGRS2280 

90  CONTINDB  D6RS2290 

GO  TO  135  DGRS2300 

95  IF  (MITER. NE.O)  GOTO  (125,125,105),  MITER  DGRS2310 

D6RS2320 

IN  THE  CASE  OF  FUNCTIONAL  ITBRATION,  06RS2330 

UPDATE  Y  DIRECTLY  FROM  THE  RESULT  DGRS2340 

OF  THE  LAST  FCTI  CALL.  D(H(S2350 

D  »  O.DO  DGRS2360 

DO  100  Ial,N  D6RS2370 

R  »  H*SAVB2(I) -Y(I,2)  DGRS2380 

D  «  D+((R-ERROR(l)  )/YMAX(l))**2  DCHIS2390 

SAVEl(I)  >  Y(I,1)+EL(1)*R  D(3RS2400 

100  ERROR  (I)  >  R  I)6RS2410 

GO  TO  145  DGRS2420 

IN  THE  CASE  OF  THE  CSIORO  METHOD,  DORS2430 

(XHPUTE  THE  CORRECnOR  ERROR,  F  SOB  OGRS2440 

(M)  ,  AND  SOLVE  THE  LINEAR  SYSTEM  DGRS2450 

WITH  THAT  AS  RIGHT-HAND  SIDS  AND  P  D6RS2460 
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AS  CX»FPICIBNT  MATRIX.  USING  THE 
LU  OBCOMPOSITIOR  IF  MITBR  ■  1  OR 
2.  IP  MITBR  -  3,  THB  COBPFICIBNT 
H*BL(1)  IN  P  IS  UPDATED. 

105  PHLO  -  HLO 
HLO  -  H*BL(1) 

IF  (HLO. BQ. PHLO)  GO  TO  115 
R  -  HLO /PHLO 
DO  110  I«1.N 

D  ■  l.D0-R*(l.D0-l.O0/PW(I) ) 

IF  (DABS(D) .BQ.O.DO)  GO  TO  165 

110  PW(I)  >  l.DO/D 

115  DO  120  I»1,N 

120  SAVEKI)  -  PW(I)*(H*SAVE2{I)-(Y(I,2)+BRRORiI))) 

GO  TO  135 

125  DO  130  I-1,N 

130  SAVEKI}  *  H*SAVE2(I)  -  (Y(I.2)-«-BRROR(I)) 

IP  (NLC  .EQ.  -1)  GO  TO  131 

NWK  *  (NLC+NUC-fl)  *N0-t-l 

CALL  LBQTlB(PH,N,NLC,NUC.N0,SAVEl,l,N0,2,PW(Minc)  ,  JER) 

GO  TO  135 

131  CALL  UJELMF  (PW, SAVBl. IPIV.N. NO, SAVED 

135  D  »  O.DO 

DO  140  I>1,M 

ERROR  (I)  -  ERROR(I)-t-SAVEl(I) 

D  «  0-t-(SAVEl(I)/YMAX(I))**2 

140  SAVEKI)  »  Y(I.l)-i-EL(l)*ERROR(I) 

TEST  FOR  CONVEROSNCE.  IF  M.GT.O,  THE 
SQUARE  OP  THB  CONVERGENCE  RATE 
CONSTANT  IS  ESTIMATED  AS  CRATE, 

AND  mis  IS  USED  IN  THE  TEST. 

145  IF  (M.NE.O)  CRATE  «  DMAXK .9D0*CRATE,D/D1) 

IF  ( (D*DMINK1.D0,2.D0*CRATE) )  .LE.BND)  GO  TO  170 
D1  «  D 
M  «  M+1 

IF  (M.EQ.3)  GO  TO  150 

CALL  FCN  (N,T,SAVEl,SAVE2,epriine,«prioie2) 

GO  TO  95 

THE  CORRECTOR  ITERATION  FAILED  TO 
CONVERGE  IN  3  TRIES.  IF  PARTIALS 
ARE  INVOLVED  BUT  ARB  NOT  UP  TO 
DATE,  THEY  ARE  REEVALUATED  FOR  THE 
NEXT  TRY.  OTHERWISE  THE  Y  ARRAY  IS 
RETRACTED  TO  ITS  VALUES  BEFORE 
PRSDICTICm,  AND  H  IS  REDUCED,  IF 
POSSIBLE.  IF  NOT,  A  NO -CONVERGENCE 
EXIT  IS  TAKEN. 

150  NFE  -  NFB+2 

IF  (IWEVAL.EQ.-l)  GO  TO  165 

155  T  «  TOLD 

RMAX  «  2. DO 
DO  160  J1>1,NQ 
DO  160  J2«J1,NQ 
J  «  (NQ+Jl) -J2 
DO  160  I>1,N 

160  Y{I,J)  -  Y(I, J) -Y(I, J+1) 

IF  (DABS(H) .LB.HMIN*1.00001D0)  GO  TO  280 
RH  «  .25D0 
IRBDO  >  1 
GO  TO  35 

165  IWFVAL  «  MITER 
GO  TO  60 
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170  IP  (MITBR.NB.O)  INBVU.  -  -1 
HFB  -  NFB-t-M 
D  •  O.DO 
DO  175  I«1,R 

175  D  •  Di-(BRROR(I)/YMAX(I))**2 
IF  (D.GT.B)  GO  TO  190 


TBB  OOniiBCTOil  KAS  OXfVBRaBD.  imVAL 
18  SIT  TO  -1  IF  PARTIAL 
D1RIVATIVB8  NIRI  USO,  TO  SIOHAL 
1VAT  THY  KMT  HID  UPDATHIG  ON 
SDBSBQOINT  8TIPS.  TH  IRROR  T88T 
IS  lADI  AND  CONTROL  PASSIS  TO 
STAmilNT  190  IF  IT  FAILS. 


KFLAG  . 
IRIDO  « 
NSTBP  a 
HDSBD  a 
BQiraBD 


0 

0 

KSTBP+l 

H 

.  HQ 


AFTH  A  StXrCHSFOL  8TBP,  DPDATB  THE 
Y  ARRAY.  CQNSIDH  CHAHOIHO  H  IF 
IDOOB  -  1.  OTHRWISI  DICRIASB 
IDODB  BY  1.  IF  IDOOB*  IS  THH  1  AND 
VQ  .LT.  MAXDH,  THH  BRROR  IS 
SAVED  FOR  USB  IN  A  POSSIBLE  ORDER 
INCRIASB  ON  TH  NEXT  STEP.  IF  A 
CHAHH  IN  H  IS  CCNSIOBRID,  AN 
HCHASI  OXt  DECREASE  IN  Oia»R  BY 
OH  IS  (Y3HIDBRBD  ALSO.  A  CHANGE 
IN  H  IS  MADE  ONLY  IF  IT  IS  BY  A 
FACTOR  OF  AT  LEAST  1.1.  IF  NOT, 
IDODB  IS  SET  TO  10  TO  PREVENT 
TESTING  FOR  THAT  MANY  STEPS. 


DO  180  Jail,L 
DO  180  Iail,N 

180  Y(I,J)  ai  Y(I,J)4BL(J)*BRRCHI(I) 
IF  (IDODB. EQ.l)  GO  TO  200 
IDODB  «  IDODB-1 
IF  (IDODB. or.  1)  GO  TO  290 
IF  (L.BQ.LMAX)  GO  TO  290 
DO  185  I-1,R 

185  Y(I,LMAX)  -  ERROR(I) 

GO  TO  290 


TH  BRROR  TEST  FAILED.  KFLAG  KEEPS 
TRACK  OF  MDLTIPLB  FAILDHS. 

HSTHB  T  AND  TH  Y  ARRAY  TO  THIR 
PRBVIOOS  VALDES,  AND  PHPAH  TO 
TRY  IH  STEP  AGAIN.  COMPUTE  TH 
OPTIMDM  STEP  SIZE  FOR  THIS  OR  OH 
LOHR  ORDH. 


190  KFLAG  -  KFLA6-1 
T  «  TOLD 
DO  195  Jlail.HQ 
DO  195  J2-J1,NQ 
J  >  (MQ+J1)-J2 
DO  195  Iail,N 

195  Y (I, J)  »  Y(I,J) -Y(1,J+1) 

RMAX  «  2. DO 

IF  (DABS(H)  .LB.HMIN«1.00001D0)  GO  TO  270 

IF  (KFLAG. LE. -3)  GO  TO  260 

IRBDO  -  2 

PR3  «  l.D-i-20 

GO  TO  210 
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SKMtDLISS  G9  TBB  8DCCB88  OR  PAIIORB 

or  TBM  am,  irctors  rri,  pr2,  mb) 

ras  ARS  OQIffORD,  BY  RRICR  H  COOLD 
BB  DZVZI»D  AT  CBtDBR  BQ  •  1,  ORtOR 
NQ,  OR  ORDBR  NQ  -f  1,  RBSRBCTIVBLY . 
ZB  TBB  CA8B  OP  FAILDRB,  PR3  • 
1.B30  TO  AVOID  AB  (BffiBR  IBCRBASB. 
IBB  SMAXiLBST  OP  TBB8B  IS 
MRBRICmD  ABD  TBB  BBW  ORDBR 
CHt^BB  AOOBtDIBSLY.  IP  TBB  OmVR 
IS  TO  BB  IBCRBASBD,  BB  COMPDTB  CRIB 
AUDinORAL  SCAL8D  DBRIYATIVB. 


200  PR3  -  l.D+20 

IP  (L.BQ.UnX)  GO  TO  210 
D1  -  O.DO 
DO  205  I-1,B 

205  D1  -  Dl+((BRROR(Z)  •Y(I,LICAX))/Y1IAX(I))**2 
BBQ3  »  .5D0/(I.>1) 

PR3  «  ((Dl/BUP)**BBQ3)*1.4D0-t-1.4D-6 
210  BBQ2  «  .5D0/L 

PR2  -  ((D/B)**BBQ2)*1.2D0+1.2D-6 
PRl  «  l.D+20 
IP  (NQ.BQ.l)  GO  TO  220 
D  «  O.DO 
DO  215  I«1,N 

215  D  -  0+(Y(I,L)/YIIAX(I))**2 
BBQl  -  .5D0/BQ 

PRl  «  ((D/BDR)**BBQ1)*1. 300+1. 3D-£ 

220  IP  (PR2.LB.PR3)  GO  TO  225 
IP  (PR3.lt. PRl)  GO  TO  235 
00  TO  230 

225  IP  (PR2.QT.PR1)  GO  TO  230 
BBBQ  •  BQ 
RH  «  1.D0/PR2 
00  TO  250 
230  BBBQ  -  BQ-1 
RH  «  l.DO/PRl 

IP  (KPLAG.BB.O.ABD.RH.GT.l.DO)  RH  «  l.DO 
00  TO  250 
235  BBBQ  «  L 

RH  ■  1.00/PR3 

IF  (RH.LT.l.lDO)  GO  TO  245 
DO  240  I>1,N 

240  Y(I,BEBQ+1)  «  BRROR(I) *BL(L) /L 
GO  TO  255 
245  IDOOB  «  10 
GO  TO  290 

250  IP  ((KFLAG.BQ.O) .ABD. (RH.LT.l.lDO))  GOTO  245 


IF  (BBBQ.BQ.NQ)  GO  TO  35 
255  NQ  «  BBBQ 
L  •  NQ+1 
IRBT  •  2 
GO  TO  15 


IF  THERE  IS  A  CHABOB  OF  ORDBR,  RESET 
BQ,  L,  ABD  THE  COBFFICIBBTS .  IB 
ABY  CASE  H  IS  RESET  ACCORDING  TO 
RH  AND  THE  Y  ARRAY  IS  RESCALED. 
TBBB  EXIT  FROM  285  IF  THE  STEP  NAS 
OK,  OR  REDO  THE  STEP  OTHERBISE. 


CONTROL  REACHES  THIS  SECTION  IF  3  OR 
MORE  FAILURES  HAVE  OCCORBD.  IT  IS 
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'  i  '  '  ' 

^  . 


ASaOBD  THAT  THE  DSRIVATIVBS  THAT 
HIVI  ACCOMDLATID  »  TOi  Y  ASKAY 
BAVi  niuRS  OF  THB  mona  exam. 
miKM  THB  FXS8T  DBSIVATIVB  IS 
RTOOMFORD,  AMD  THB  ORDER  IS  SBT 
TO  1.  TRBH  K  IS  RBOOCBD  BY  A 
FACTCR  OF  10,  AMD  THB  STEP  IS 
RETRIED.  AFTER  A  TOTAL  OF  7 
FAILraOtS,  AH  EXIT  IS  TAXBlf  WITH 
KFLAO  -  -2. 

260  IF  (KFLAO. BQ. -7)  (30  TO  275 
mt  m  .IDO 

im  >  DMAXl  (HMIN/DABS  (H) ,  RH) 

H  ■  H*RH 

CALL  FOI  (B,T,Y.SAVBl,«prin«,«priaM2) 

HFE  -  RFE41 
DO  265  I«1,R 
265  Y(I,2)  -  H*SAVB1(I) 

IHBVAL  *  mTBR 
IDOOB  -  10 

IF  (NQ.BQ.l)  00  TO  SO 
HQ  •  1 
L  •  2 
IRBT  «  3 
GO  TO  IS 

ALL  RETORHS  ARE  MADE  THROIXS  THIS 
SBCTIQH.  H  IS  SAVED  IN  BOLD  TO 
ALLOW  THE  CALLER  TO  CBAMGB  H  ON 
THB  NEXT  STEP. 

270  KFLAO  •  -1 
GO  TO  290 
275  KFLAO  •  -2 
GO  TO  290 
280  KFLAO  -  -3 
GO  TO  290 
285  RMRX  •  10. DO 
290  HOLD  «  H 

JSTART  -  HQ 

C* -Diagnostic  Chock  of  first  and  sacond  derivativaa  of  B 
if (tcum.aq.told)9o  to  310 

writs (8, 300) tcum, atop, y (1,1) ,y(2,l) ,y(3,l) ,y(4,l) ,y(5,l) 

300  foniiat(lx,ell.4,lx,x5,5(lx,ell.4)) 
write  (9,305)  step,  eprine ,  epriiae2 
305  fonoat (lx, 15,2 (Ix,e20.l3) ) 

RETORH 

END 
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mSL  SOOTZIIB  NMOI 


IX»CS 


coMErariR 
LKIBST  RKVISICW 
FORPOSB 

PSBCISiai/IiARDIlAltB 

RBQD.  XUSL  ROtmilBS 
NOTftnOR 


•  IBH/DOOBUS 

•  JANOARy  1,  1978 

•  NOCLBOS  CALLED  (»ILY  BY  IMSL  SDBROOTINE  DGBAR 

•  SmQLB  ABD  D0DBLB/H32 

•  SIMGLB/Ha6,H48,H€0 


COPYRIGHT 

KARRANTY 


MORE  RBgOIRSD 

•  INFORUATKm  ON  SPECIAL  NOTATICM  AND 

COHVENTICRIS  IS  AVAILABLE  IN  THE  MANOAL 
INTRODOCnON  OR  THROOQH  IMSL  RODTINE  UHELP 

•  1978  BY  IMSL,  INC.  ALL  RKSITS  RESERVED. 

-  IMSL  WARRANTS  ONLY  THAT  IMSL  TESTING  HAS  BEEN 
APPLIED  TO  THIS  CODE.  NO  OTHER  WARRANTY, 
EXPRESSED  OR  IMPLIED,  IS  APPLICABLE. 


SOBROOTINE  DGRCS 


INTEGER 

REAL 

DODBLE  PRECISICOI 

INTEGER 

REAL 

DATA 


(MBTH,  NQ,  EL,TQ,MAXDER) 

SPECIPICATIONS  FOR  ARGOMENTS 
METH,NQ,MAXDER 
TQ(1) 

EL(1) 

SPECIFICATIONS  FOR  LOCAL  VARIABLES 
K 

PERTST(12,2,3) 

PEHTST/1.,1.,2.,1., .3158, .7407B-1, 

1  .1391B-1, .2182B-2, .2945B-3, .3492B-4, 

2  .3692B-5, .3524B-6,1. ,1. , .5, .1667, 

3  . 4167B- 1 , 7*1 , , 2 . , 12 . , 24 . , 37 . 89 , 

4  53.33,70.08,87.97,106.9,126.7, 

5  147.4,168.8,191.0,2.0,4.5,7.333, 

6  10. 42, 13. 7, 7*1., 12. 0,24. 0,37. 89, 

7  53.33,70.08,87.97,106.9,126.7, 

8  147. 4, 168. 8, 191. 0,1. ,3. 0,6.0, 

9  9.167,12.5,8*1./ 

FIRST  BXECDTABLB  STATEMENT 

GO  TO  (5,10) ,  METH 
5  MAXDER  «  12 

GOTO  (15,20,25,30,35,40,45,50,55,60,65,70),  NQ 
10  MAXDER  .  5 


GOTO  (75,80,85,90,95),  NQ 


THE  FOLLOWING  COEFFICIENTS  SHOULD  BE 
DEFINED  TO  MACHINE  ACCURACY.  FOR  A 
GIVEN  ORDER  NQ,  THEY  CAN  BE 
CALCULATED  BY  USB  OF  THE 
GENERATING  POLYNOMIAL  L  (T)  ,  WHOSE 
COEFFICIENTS  ARB  EL(I) . .  L(T)  « 
EL(1)  +  BL(2)*T  +  . . .  + 
EL(NQ+1)*T**NQ.  FOR  THE  IMPLICIT 
ADAMS  METHODS,  L(T)  IS  GIVEN  BY 
DL/DT  .  (T+1) * (T+2) *  . . . 
*(T+NQ-1)/K,  L(-l)  -  0,  WHERE  K  = 
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c 

c 

c 

c 

c 

c 

c 

c 

c 


MCTC»IAL(aQ-l)  .  POR  tHB  OSAR 
IBTBOOS,  L(T)  -  (T«l)*(T-f2)*  ... 

*(T4SQ)/k,  whirs  K  - 

VACTORZAX.(MQ)*(l  *  X/2  *  ... 
l/MQ)  .  THS  (XnSR  IS  WIICH  TRB 
ORODRS  APPSAR  BSLON  IS.  .  IMPLICIT 
A&AMS  1BTHQ08  OP  ORZORS  1  TO  12, 
BACKKARD  DIPFSRBMriATIOM  MSTMODS 
OF  ORDSRS  1  TO  5. 

15  BL(1)  -  1.000 
GO  TO  100 

20  SL(1)  •  0.500 
SL(3)  -  0.500 
GO  TO  100 

25  SL(1)  •  4.1S€66666e666667D-01 
SL(3)  -  0.7500 

BL(4)  ■  1.6S€666S£66€56fi70-01 
GO  TO  100 

30  BL(1)  >  0.37500 

BL(3)  «  9.16£€€£6£€£S6€€70-01 
BL(4)  •  3.3333333333333330-01 
BL(5)  -  4.1£££££££££££££70-02 
GO  TO  100 

35  BL(1)  «  3.48fillllllllllllO-01 
SL(3)  «  1.041£££££££££££700 
BL(4)  >  4.8£11111111111110-01 
BL(5)  -  1.041£££££££££££70-01 
BL(£)  •  8.3333333333333330-03 
GO  TO  100 

40  BL(1)  «  3.298£111111111110-01 
SL(3)  -  1.141£££££££££££7D-«-00 
BL(4)  -  0.8250+00 
SL(S)  -  1.7708333333333330-01 
SL(£)  «  0.0250+00 
SL(7)  .  1.3888888888888890-03 
GO  TO  100 

45  BL(1)  «  3.1S5919312189312D-01 
BL(3)  -  1.2250+00 
Bu(4)  «  7.5185185185185190-01 
BL(5)  -  2.5520833333333330-01 
HL(6}  -  4.881111111111111O-02 
BL(7)  «  4.881111111111111O-03 
BL(8)  *  1.98412€98412£984D-04 
GO  TO  100 

50  BL(1)  -  3.0422453703703700-01 
BL(3)  >  1.2984285714285710+00 
BL(4)  *  8.8851851851851850-01 
BL(5)  »  3.3578388888888890-01 
BL(8)  -  7.7777777777777780-02 
BL(7)  «  1.0848148148148150-02 
BL(8)  -  7.9385079385079370-04 
BL(9)  -  2.4801587301587300-05 
GO  TO  100 

55  BL(1)  «  2.9488800044091710-01 
BL(3)  s  1.3589285714285710+00 
BL(4)  »  9.7855423280423280-01 
BL(5)  >  4.1718750-01 
BL(6}  «  1.1135418866686670-01 
BL(7)  «  0.018750+00 
BL(8)  -  1.9345238095238100-03 
BL(9)  -  1.1160714285714290-04 
BL(10)b  2.7557319223985890-06 


0(BtC0830 

OQRC0840 

»»C0850 

PGRC0880 

3abC0870 

OQRC0880 

OaiC0890 

0(MtC0700 

0(MU:0710 

OCBtC0720 

OGRC0730 

OGntC0740 

06RC0750 

0(MtC0760 

DffitC0770 

OGRC0780 

OGRC0790 

OGRC0800 

OGRC0810 

OGffiC0820 

0(MtC0830 

0«IC0840 

OGRC0850 

OGRC0860 

0^0870 

OCBtC0880 

omcoaso 

OGRC0900 

OGRC0910 

OGRC0920 

00100930 

OGRC0940 

OGRC0950 

00100980 

OGRC0970 

00100980 

OGRC0990 

OORCIOOO 

OGRClOlO 

OQRC1020 

OGRC1030 

OGRC1040 

00«C1050 

OGRC1080 

OGRC1070 

OGRC1080 

06RC1090 

DGRCllOO 

OGRClllO 

OGRC1120 

OGRC1130 

OGRC1140 

OGRC1150 

OORC1180 

OGRC1170 

OGRC1180 

OGRC1190 

06RC1200 

OGRC1210 

OGRC1220 

OGRC1230 

OGRC1240 
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00  lO 

100 

OOItC1350 

60  IL(1) 

2.8697544642857140-01 

DatC1260 

BX.(3) 

1.4144841269841270+00 

O(»C1270 

0.(4) 

1 . 0772156084656090+00 

DGItC1280 

BL(5) 

4.9856701940035270-01 

D6RC1290 

IL(6) 

1.4843750-01 

DaRC1300 

IL(7) 

2.9060570987654320-02 

OGSC1310 

BL(8) 

3.7202380952380950-03 

DGRC1320 

IL(9) 

2.9968584656084660-04 

D(aC1330 

IL(IO) 

1 . 3778659611992950- 05 

DGmC1340 

BL(IX) 

2.7557319223985890-07 

DOItC1350 

GO  to 

100 

DORC1360 

65  BL(1) 

2 . 801895964439367D-01 

DGIU:i370 

BL(3) 

1.4644841269841270+00 

DGRC1380 

BL(4) 

1 . 1715145502645500+00 

DaRC1390 

BL(5) 

5.7935819003527340-01 

DGRC1400 

BL(6) 

1 . 8832286155202820-01 

OORC1410 

BL(7) 

4 . 1430362654320990-02 

DGRC1420 

EL(8) 

6.2111441798941800-03 

DGRC1430 

BL(9) 

6.2520667989417990-04 

DGRC1440 

BL(IO) 

4.0417401528512640-05 

DGRC14S0 

BLdl) 

1.5156525573192240-06 

DGRC1460 

BL(12) 

2 . 5052108385441720-08 

DGRC1470 

GO  TO 

100 

OGRC1480 

70  BL(1) 

m 

2 . 7426554003159910-01 

DGRC1490 

BX.(3) 

m 

1 . 5099386724386720+00 

DORC1500 

BL(4) 

m 

1 . 2602711640211640+00 

DORC1510 

BL(5) 

m 

6.5923418209876540-01 

D(»IC1520 

BL(6) 

m 

2.3045800264550270-01 

DORC1530 

BL(7) 

m 

5 . 5697246105232220- 02 

DGRC1S40 

B2<(8) 

m 

9.4394841269841270-03 

DGRC1S50 

BL(9) 

m 

1.1192749669312170-03 

DGRC1560 

BL(IO) 

«  9.09391S343915344D-05 

DGRC1570 

BL(ll) 

>  4.822530864197S31D’06 

OGatC1580 

BL(12) 

m 

1.5031265031265030-07 

DGRC1590 

BL(13) 

m 

2.0876756987868100-09 

DGRC1600 

GO  TO 

100 

DGRC1610 

DGRC1620 

75  BL(1) 

m 

1.00+00 

DGRC1630 

GO  TO 

100 

DGRC1640 

80  BL(1) 

m 

6.66666666€  670-01 

DGRC1650 

BL(3) 

M 

3.33333333  53  330-01 

DGRC1660 

GO  TO 

100 

DGRC1670 

85  BL(1) 

m 

5.4545454545454550-01 

DGRC1680 

BL(3) 

m 

BLd) 

DGRC1690 

BL(4) 

m 

9.0909090909090910-02 

OGRC1700 

GO  TO 

100 

DGRC1710 

90  BX.(1) 

s 

0.480+00 

DGRC1720 

BL(3) 

■ 

0 . 70+00 

DGRC1730 

BL(4) 

m 

0.20+00 

OGRC1740 

BL(5) 

0.020+00 

DGRC1750 

GO  TO 

100 

DGRC1760 

95  BL(1) 

m 

4.3795620437956200-01 

DGRC1770 

BL(3) 

m 

8.2116788321167880-01 

DGRC1780 

BL(4) 

m 

3.1021897810218980-01 

DGRC1790 

BL(5) 

a 

5.4744525547445260-02 

DGRC1800 

BI.(6) 

a 

3 . 6496350364963500-03 

DGRC1810 

DGRC1820 

100  DO  105  K«l,3 

DGRC1830 

TQ(K)  -  PBRTST(NQ,METH,K) 

DGRC1840 

105  CONTINDB 

DGRC1850 

TQ(4) 

a 

.5D0*TQ(2)/(HQ+2) 

DGRC1860 
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DOICISTO 

DCRC1880 
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ZMSL  ROOTIHB 


DORPS 


COMDPCITBR 
XATBST  RBVISIC» 
PORPOSB 

PRBCISZQR/HARDWARB 


-  IBM/DOOBLB 

•  NOVBMBBR  1,  1984 

•  NOCLBOS  CAT.T.Bn  ORLY  BY  IMSL  SUBROCTZHB  DGBAR 

-  SZNQLB  AND  DOaBLB/H32 

•  SINOLB/H36,H48,H60 


RBQD.  ZMSL  RODTZNBS  -  LDDATF,LBQTlB,DBRTST,OGBTZO 

NOTATZCBH  •  ZRFORMATZON  ON  SPBCZAL  NOTATZON  AND 

CONVBNTZONS  ZS  AVAZLABLB  ZN  THB  MAHOAL 
ZNTRODOCTZON  OR  THROOGOl  ZMSL  ROOTZNB  DHBLP 

CX)PYRZGHT  -  1984  BY  ZMSL,  ZNC.  ALL  RZGHTS  RBSBRVBD. 

WARRANTY  -  ZMSL  WARRANTS  ORLY  THAT  ZMSL  TBSTZNG  HAS  BBBN 

APPLZBD  TO  THZS  C(B}B.  NO  OTHBR  WARRANTY, 
BXPRBSSBD  OR  ZMPLZBD,  ZS  APPLZCABLB. 


SDBROOTZNB  DORPS 


ZNTBOBR 

DODBLB  PRBCZSZON 


ZNTBOBR 


REAL 

DODBLB  PRBCZSZCXI 
COMMON  /DBAHD/ 

caiaaost  /obar/ 


(FCN,  FCRJ,  Y,  NO ,  CON,  MZTBR,  YMAX,  SAVBl ,  SAVB2 ,  PW, 
BQDZL.ZPZV,  ZBR) 

SPBCZFZCATZORS  FOR  ARaOMBNTS 
NO, MZTBR,  ZPZV(l)  ,ZBR 

Y  (NO,  l),C(»r,  YMAX  (1),  SAVBl  (1)  ,SAVB2  (1)  ,PW(1) , 
BQOZLd) 

SPBCZFZCATZONS  FOR  LOCAL  VARZABLBS 

NC ,  MFC , KFLAG , JSTART , NQDSBD , NSTBP , NFB , NJB , NPW , 
HSQ, Z , J1 , J, NBRROR, NSAVBl , RSAVB2 , NBQDZL, NY, 
ZDDMIY(23)  ,NLZM,ZZ,ZJ,LZMl,LZM2,RB,NLC,NOC,NWK 
SDDIMY(4) 

T, H, HMZN, HMAX, BPSC, DRODRD , BPS J , HDSBD , D , RO , YJ , R, 
D1 ,  D2 ,  WA,  DDMMY  (40 ) 

NLC,NDC 

T , H , HMZN , HMAX , BPSC , URODND , BPS J , HDSBD , DUMMY , 
SDDMCY ,  NC ,  MFC ,  KFLAG ,  JSTART ,  NSQ ,  NQDSBD ,  NSTBP , 
NFB , NJB , NPW, NBRROR, NSAVBl , NSAVE2 , NBQDZL, NY, 
ZDDMIY 

THZS  RODTZNB  ZS  CALLBD  BY  DGRST  TO 
COMPDTB  AND  PROCBSS  THB  MATRZX  P  > 
Z  -  H*BL(1)*J  ,  WHBRE  J  ZS  AN 
APPROXZMATZON  TO  THB  JACOBZAN.  J 
ZS  COMPDTBD,  EZTHER  BY  THE  DSER- 
SDPPLZED  RCXZrZNE  FCNJ  ZF  MZTBR  - 
1,  OR  BY  FZNZTB  DZFFBRBNCZNG  ZF 
MZTBR  .  2 .  J  ZS  STORED  ZN  PW  AMD 
REPLACED  BY  P,  DSZMG  CON  > 

•H*BL(1) .  THEN  P  ZS  SUBJECTED  TO 
LD  DECXUPOSZTZOM  ZN  PREPARATION 
FOR  LATBR  SOLDTZON  OF  LZNBAR 
SYSTEMS  WZTH  P  AS  COBFFZCZENT 
MATRZX.  ZN  ADDZTZOM  TO  VARIABLES 
DBSCRZBED  PREVZODSLY, 

COMMDMZCATZOM  WZTH  DGRPS  USES  THE 


DQtPOOlO 

DOtPOOZO 

-D«P0030 

DORP0040 

D(»P0050 

DCBU>0080 

DORP0070 

DORP0080 

D6RP0090 

DORPOlOO 

DORPOllO 

DCHtPOlZO 

DGRP0130 

D(MU>0140 

DCRPOISO 

DORP0160 

DORP0170 

DGRP0180 

D(3RP0190 

DORP0200 

DORP0210 

DGRP0220 

OGRP0230 

DGRP0240 

DORP0250 

-DORP0260 

DORP0270 

DCRP0280 

DORP0290 

DOIPOZOO 

DGRP0310 

DORP0320 

DORP0330 

DCHIP0340 

DORP0350 

DGRP0360 

OORP0370 

DGRP0380 

DGRP0390 

DORP0400 

DGRP0410 

DGRP0420 

DORP0430 

DORP0440 

DGRP0450 

DGRP0460 

DGRP0470 

DGRP0480 

DGRP0490 

DGRP0500 

DGRP0510 

DGRP0520 

DGRP0530 

DGRP0540 

DGRP0550 

06RP0560 

DGRP0570 

DGRP0580 

DGRP0590 

DORP0600 

DGRP0610 

DGRP0620 


71 


no  n  non  00000 


C 

C 


C 


roLumsa  rpsj  -  dsqrt  (orodnd)  , 

USED  nt  THB  NDMIRICAL  JACQBIAH 
ZlKntBiaHTS. 


FIRST  BXBCDTABLB  STATBMBNT 

IF  (NLC.BQ.-l)  GO  TO  45 

BANDSD  JACOBIAN  CA5B 

NB  -  NLC-t-NDC-fl 

NVK  >  NB*N0+1 

IF  (MITBR.BQ.2)  GO  TO  15 

MITBR  .  1 

NLIM  -  NB*N0 
DO  5  I-1,NLIM 
PVd)  «  O.ODO 
5  COHTIHUB 

CALL  FCNJ(NC.T,Y,I>N) 

DO  10  I>1.NLIM 

PW(I)  -  PW(I)*CON 
10  CONTINOB 
GO  TO  35 

MITBR  .  2 

15  D  ■  0.000 
DO  20  I«1,NC 
20  D  «  D-fSAVB2(I)**2 

RO  -  DABS(H)*DSQRT(D)*1.0D-«-03*URODRD 
DO  30  J-1,MC 
YJ  «  Y(J,1) 

R  *  BPSJ*YMAX(J) 

R  «  raiAZKR.RO) 

Y(J,1)  «  Y(J,1)+R 

O  •  con/R 

CALL  FCN(MC,T,Y,SAVB1) 

Lna  «  MRXOd.  j-Hoc) 

LIM2  «  MIN0(N0,  J+NLC) 

DO  25  I>LIK1,LIM2 

IJ  -  (J-I+NLC)*NO-i-I 
PNdJ)  «  (SAVBld)  •SAVB2d))*D 
25  CONTINOB 

Y(J,1)  -  YJ 
30  CONTINOB 

ADD  IDENTITY  MATRIX. 

35  DO  40  I>1,NC 

II  -  NLC*N0-t-I 
PWCII)  «  PWdD-t-l.ODO 
40  CONTINOB 

DO  m  DBCOMPOSiTicnr  on  p 

CALL  LBQTIB  (PW, NC, NLC,  NOC, NO , EQOIL,  1,N0 , 1 ,  PW (NWK)  ,  lER) 
RBTORN 

POLL  JACOBIAN  CASE 

45  IF  (MITBR. EQ. 2)  GO  TO  55 

MITER  .  1 

CALL  FCNJ(NC,T,Y,PW) 

DO  50  I>1,NSQ 
50  PWd)  «  PWd)*CON 
GO  TO  75 

MITBR  «  2 

55  D  -  O.ODO 
DO  £0  I>1,NC 
60  D  «  Di-SAVE2d)**2 

RO  >  DABS(H)*DSQRT(D)*1.0D4-03«ORODND 
J1  >  0 


DW0630 

DQRP0640 

DCnP0650 

DCBIP0660 

D6RP0670 

DORPOSSO 

OaRP0690 

DGRP0700 

DGRP0710 

D(»P0720 

DaRP0730 

D6RP0740 

DGRP0750 

DGRP0760 

DaRP0770 

DCAP0780 

DGniP0790 

DGRPOSOO 

DGRP0810 

OGRP0820 

DaRP0830 

DGRP0840 

DGRP0850 

DGRP0860 

DORP0870 

DQRP0880 

D(»P0890 

DORP0900 

D(atP0910 

D(»P0920 

DaRP0930 

DGRP0940 

OGRP0950 

DORP0960 

DGRP0970 

0au>0980 

DaRP0990 

DGRPIOOO 

DGRPlOlO 

DratP1020 

DORP1030 

DGRP1040 

DGRP1050 

DORP1060 

D6RP1070 

DGRP1080 

DGRP1090 

DGRPllOO 

DGRPlllO 

DGRP1120 

DORP1130 

DGRP1140 

DGRP1150 

DGRPH60 

D^1170 

D^llSO 

DGRP1190 

DORP1200 

D6RP1210 

DGRP1220 

DGRP1230 

D6RP1240 
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DO  70  J«1,NC 
YJ  .  Y(J,X) 

R  >  BPSJ*YMAX(J) 

R  -  DMAXKR.RO) 

Y(J,1)  -  Y(J,1)+R 
D  «  ccai/R 

CALL  FCN(NC.T,Y,SAVB1) 

DO  65  I-l.NC 

65  PW(I-t-Jl)  «  (SAVBl(Z) -SAVB2(I))«D 
Y(J,1)  -  YJ 
J1  -  Jl+NO 
70  CONTINOB 

C  ADD  IDENTITY  MATRIX. 

75  J  -  1 

DO  80  I>1,NC 

PW(J)  «  PW(J)-t-1.0D0 
J  s  J+(N0-I>1) 

80  CONTINUE 

C  DO  LU  DECOMPOSITION  ON  P. 

C 

CALL  LUDATF (PW, PW, NC , NO , 0 . D1 , D2 , IPIV, BQDIL, HA, lER) 

RETURN 

END 


DQRP1250 

DORP1260 

D6RP1270 

DGRP1280 

DGRP1290 

DGRP1300 

DQRPiaiO 

DGRP1320 

DGRP1330 

DGRP1340 

DQRP1350 

DQRP1360 

D6RP1370 

DGRP1380 

DQRP1390 

DGRP1400 

D6RP1410 

D6RP1420 

DGRP1430 

D6RP1440 

DGRP1450 

DGRP1460 

DGRP1470 
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onnooofjnonn  o  no  oonoonooooononoooonnooonoon 


niSL  RODTINB  NAMB  -  DGRIN 


COMFOTHR 
LKTBST  RBVISICKf 
ramosB 

PRBCISICm/HARDHARB 


IBM/DOCBLB 
JANOARY  1.  1978 


.  NDCLEOS  CALLBD  ONLY  BY  IMSL  SUBROOTINE  DGBAR 

•  SINGLE  AND  DOaBLB/H32 
-  SINQLE/H36,K48,H60 


REQD.  IMSL  ROUTINES  -  NONE  REQUIRBO 


NOTATION 


CX3PYRI6HT 


NARRANTY 


INFORMATION  ON  SPECIAL  IKITATION  AND 
CONVENTICBSS  IS  AVAILABLE  IN  THE  MANUAL 
INTRCn>nCTION  OR  THRCXm  IMSL  ROUTINE  UHBLP 

1978  BY  IMSL,  INC.  ALL  RIGHTS  RESERVED. 


DQRIOOlO 

DGRI0020 

•Dcgtiooao 

D(nu[0040 

DQRI0050 

DGRI0060 

D(Ma0070 

DGRIOOSO 

DGRI0090 

D6RI0100 

DGRIOllO 

DGRI0120 

D6RI0130 

DGRI0140 

DGRI0150 

DGRI0160 

DQRI0170 

DGRI0180 

D6RI0190 

OGRI0200 


SUBROUTINE  DGRIN 
INTEGER 

DOUBLE  PRECISION 
INTEGER 


REAL 

DOUBLE  PRECISION 
L 

COMMON  /GEAR/ 


DO  5  I  »  1,NC 

YO(I)  =  Y(I,1) 
5  CONTINUE 


DORI0210 

•  IMSL  WARRANTS  ONLY  THAT  IMSL  TESTING  HAS  BEEN  DGRI0220 

APPLIED  TO  THIS  CODE.  NO  OTHER  WARRANTY,  D6RI0230 

EXPRESSED  OR  IMPLIED,  IS  APPLICABLE.  DGRI0240 

DGRI0250 

. D(»tI0260 

imi0270 

(TOUT,Y,NO,YO)  DQRZ0280 

SPECIFICATIONS  FOR  ARGUMENTS  DCnRI0290 

NO  DGRI0300 

T0IJT,Y0(N0)  ,Y(N0,1)  DGRI0310 

SPECIFICATKWS  FOR  LOCAL  VARIABLES  D6RI0320 

NC,MFC,KFLA6,I,L,  J,  JSTART,NSQ,NQaSBD,NSTBP,  OGRI0330 

NFB,NJE,NPW,NBRR0R,NSAVB1,NSAVB2,NBC2UIL,NY,  DGRI0340 

IDUMMY(23)  DCBII0350 

SO0UMY(4)  D6RI0360 

T,H,HMIN,HMAX,BPSC, GROUND, BPSJ,HUSBD,S, SI,  DORI0370 

DUMMY (40)  DGRI0380 

T,H,HMIN,HMAX,BPSC,UROUMD,BPSJ,HUSED,DUiafY,  DGRI0390 

SDUiaiY,NC,MFC,KFLA6,JSTART,NSQ,NgUSED,NSTBP,  DGRI0400 

NFE , NJE , NPW, NERROR, NSAVEl , NSAVE2 , NEQUIL, NY,  DGRI0410 

IDUMHY  DGRI0420 

FIRST  EXECUTABLE  STATEMENT  OGRI0430 

DGRI0440 
DGRI0450 
OGRI0480 

THIS  SUBROUTINE  CCBIPUTES  INTERPOLATEDDGRI0470 


L  >  JSTART  +  1 
S  «  (TOUT  -  T)/H 
SI  «  l.ODO 
DO  15  J  -  2, L 
'SI  >  S1*S 


VALUES  OF  THE  DEPENDENT  VARIABLE 
Y  AND  STORES  THEM  IN  YO.  THE 
INTERPOLATION  IS  TO  THE 
POINT  T  s  TOUT,  AND  USES  THE 
NORDSIECX  HISTORY  ARRAY  Y,  AS 
FOLLOWS. . 

NQ 

Y0(I)  .  SUM  Y(I, J+1)*S**J  , 

J>0 

WHERE  S  =  -(T-TOGT)/H. 


DGRI0480 

OGRI0490 

DGRI0500 

DGRI0510 

DGRI0520 

DGRI0530 

DGRI0S40 

DGRI0550 

DGRI0560 

DORI0570 

DGRI0S80 

OGRI0S90 

DGRI0600 

DGRI0610 

DGRI0620 
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DO  10  I  «  1,I1C  IX3R10630 

Y0(I)  «  Y0(X)  +  S1*Y(Z,J)  DQRI0640 

10  CONTINUB  D(SU0€50 

15  CQMTIHUB  D(»II0660 

RTTOKN  DGR10670 

BBS)  D6RI0680 


75 


nnnnooooonnoooooooooonoooonooonoonnnooooonnonnooonnoooooooooon 


umKrr 


COMWTrgR 
LATIST  RSVISIOR 
PORPOSB 


-  m/DOOBIOl 

•  JMnavRy  i,  i978 

•  L-n  DBCGNPOSITION  BY  THB  CROtTr  AL60RIT»C 
WITH  OPTIOKBL  ACCURACY  TBST. 


USAGE 


CALL  LUDATP  (A.LU,N,lA,IDaT.Dl,D2,IFVT, 
BQDIL.KA.IBR) 


ARGUMENTS 


A  •  INPUT  MATRIX  OP  DIMBR8I0H  N  BY  N  C(»ITAINING 

THE  MATRIX  TO  BE  DECOMPOSED. 

LU  *  REAL  OUTPUT  MATRIX  OP  DIMENSION  N  BY  N 
CQNTAINIMO  THE  L*U  DBCOMPOSITICNI  OP  A 
ROWWISE  PERMDTATIOM  OP  THE  INPUT  MATRIX. 
POR  A  DESCRIPTION  OP  THB  PORMAT  OP  LU,  SEE 
EXAMPLE. 

N  -  INPUT  SCAxjul  CONTA1NIBK3  THE  ORDER  OP  THB 

MATRIX  A. 

lA  •  INPUT  SCALAR  CXXrEAININQ  THB  RON  DIMENSION  OP 
MATRICES  A  AND  LU  EXACTLY  AS  SPBCIPIBD  IN 
THB  CALLING  PROGRAM. 

IDGT  -  INPUT  OFTICRI. 

IP  IDGT  IS  GREATER  THAN  ZERO,  THE  WXI-ZERO 
ELEMENTS  OP  A  ARE  ASSUMED  TO  BE  CORRECT  TO 
IDGT  DECIMAL  PLACES.  LODATP  PERFORMS  AN 
ACCURACY  TBST  TO  DETERMINB  IP  THB  COMPUTED 
DECOMPOSITION  IS  THB  EXACT  DECOMPOSITION 
OP  A  MATRIX  WHICH  DIPPERS  PROM  THB  GIVEN 
ONE  BY  LESS  THAN  ITS  UNCERTAINTY. 

IP  IDGT  IS  EQUAL  TO  ZERO,  THE  ACCURACY  TBST 
IS  BYPASSED . 

D1  -  OUTPUT  SCALAR  CONTAINING  ONE  OP  THE  TWO 
COMPONENTS  OP  THB  DETERMINART.  SEE 
DESCRIPTION  OP  PARAMETER  D2,  BELOW. 

D2  •  OUTPUT  SCALAR  CONTAINZNG  ONE  OP  THB 

TWO  COMPONENTS  OP  THB  DETERMINANT.  THB 
DETERMIRANT  MAY  BE  EVALUATED  AS  (Dl)  (2**D2} 
IPVT  •  OUTPUT  VECTOR  OP  LENGTH  N  CCEYTAINING  THE 
PERMUTATION  INDICES.  SEE  DOCUMENT 
(ALGORITHM)  . 

BQUIL  -  OUTPUT  VECTOR  OP  LENGTH  N  CONTAINING 

RECIPROCALS  OF  THE  ABSOLUTE  VALUES  OP 
THB  LARGEST  (IN  ABSOLUTE  VALUE)  ELEMENT 
IN  BACH  RON. 

WA  -  ACCURACY  TEST  PARAMETER,  OUTPUT  ORLY  IF 
IDGT  IS  GREATER  THAN  ZERO. 

SEE  ELEMENT  DOCUMENTATION  FOR  DETAILS. 

IBR  •  ERROR  PARAMETER.  (OUTPUT) 

TERMINAL  ERROR 

lER  -  129  INDICATES  THAT  MATRIX  A  IS 
ALGORITHMICALLY  SINGULAR.  (SEE  THB 
CHAPTER  L  PRELUDE) . 

WARNING  ERROR 

IBR  «  34  INDICATES  THAT  TFE  ACCURACY  TBST 
FAILED.  THE  COMPUTED  SOLUTION  MAY  BE  IN 
ERROR  BY  MORE  THAN  CAN  BE  ACCOUNTED  POR 
BY  THB  UNCERTAINTY  OF  THE  DATA.  THIS 


LODAOOlO 

LDDA0020 

•LUDA0030 

.LUDAD040 

UUDAOOSO 

UE3A0060 

LDDA0070 

UE3A0080 

LODA0090 

LODAOlOO 

LUDAOllO 

LDDA0120 

LDDA0130 

LDDA0140 

LODA0150 

L0DA0160 

LODA0170 

LDDA0180 

LUDAOISO 

LODA0200 

L0DA0210 

LUDA0220 

IJ0DA0230 

LDDA0240 

LODA0250 

LODA0260 

LDDA0270 

LODA0280 

L0DA0290 

LDDA0300 

LDDA0310 

LDDA0320 

LODA0330 

L0DA0340 

LDDA0350 

LDDA0360 

LDDA0370 

LDDA0380 

LDDA0390 

LODA0400 

LODA0410 

.LUDA0420 

L0DA0430 

LDDA0440 

L0DA0450 

LnDA04€0 

LDDA0470 

IJDDA0480 

LODA0490 

LDDA0500 

LUDAOSIO 

LODA0520 

LUDA0530 

LUDA0540 

LODA0550 

LUDA0560 

LODA0570 

LODA0580 

L0DA0590 

LUDA0600 

LUDAOSIO 

LDDA0620 
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no  n  o o o o n n o n n n o n n n o o o n o o o o o n o 


NARNZMG  CAN  BE  PRODUCED  ONLY  IF  IDGT  IS 
GREATER  THAN  0  ON  INPUT.  SEE  CHAPTER  L 
PRELUDE  FOR  FURTHER  DISCUSSI(»I. 

PRECISION/HARDWARE  -  SINGLE  AND  DOOBLE/H32 

-  SINGLE/H36,H48,H60 

REQD.  I2fSL  ROUTINES  -  UERTST,UGETIO 

IR3TATION  -  INFORMATION  ON  SPECIAL  NOTATION  AND 

CONVENTIONS  IS  AVAILABLE  IN  THE  MANUAL 
INTRODUCTION  OR  THIU3UGH  IMSL  ROUTINE  UHELP 

REMARKS  A  TEST  FOR  SINGULARITY  IS  MADE  AT  TWO  LEVELS: 

1.  A  ROW  OF  THE  ORIGINAL  MATRIX  A  IS  NULL. 

2 .  A  COLUMN  BECOMES  NULL  IN  THE  FACTORIZATION  PROCESS 

COPYRIGHT  -  1978  BY  IMSL,  INC.  ALL  RIGHTS  RESERVED. 

WARRANTY  -  IMSL  WARRANTS  ONLY  THAT  IMSL  TESTING  HAS  BEEN 

APPLIED  TO  THIS  CODE.  tK)  OTHER  WARRANTY, 
EXPRESSED  OR  IMPLIED,  IS  APPLICABLE. 


SUBROUTINE  LUDATF  (A,LU,N, lA, IDGT,D1,D2, IPVT.EQUIL, WA, lER) 

DIMENSION  A(IA,1) ,LU(IA,1) ,IPVT(1) ,EQUIL(1) 

DOUBLE  PRECISION  A,LU,Dl,D2,EQUIL,WA,ZERO,ONE,FOUR,SIXTN,SIXIH, 

•  RN,WREL,BIGA,  BIG,  P,  SUM,  AI,WI,T,  TEST,  Q 

DATA  ZERO, ONE , FOUR, SIXTH, SIXTH/ 0 . DO , 1 . DO , 4 .DO , 

*  16.D0, .0625D0/ 

FIRST  EXECUTABLE  STATEMENT 
INITIALIZATION 

lER  »  0 
RN  -  N 
WREL  ZERO 
D1  «  ONE 
D2  s  ZERO 
BIGA  s  ZERO 
DO  10  lal.N 
BIG  =  ZERO 
DO  5  J=1,N 
P  -  A(I,J) 

LU(I,J)  >  P 
P  s  DABS(P) 

IF  (P  .GT.  BIG)  BIG  =  P 
5  CONTINUE 

IF  (BIG  .GT.  BIGA)  BIGA  =  BIG 
IF  (BIG  .EQ.  ZERO)  GO  TO  110 
EQUIL(I)  =  ONB/BIG 
10  CONTINUE 

DO  105  J-1,N 
JMl  s  J-1 

IF  (JMl  .LT.  1)  GO  TO  40 

COMPUTE  D( I, J) ,  1=1,..., J-1 

DO  35  1=1,  JMl 
SUM  =  LU(I,J) 

IMl  >  I-l 

IF  (IDGT  .EQ.  0)  GO  TO  25 

WITH  ACCURACY  TEST 

AI  =  DABS (SUM) 


LODA0630 

LaDA0640 

L0DA0650 

LODA0660 

LODA0670 

LDDA0880 

LUDA0690 

LUDA0700 

LDDA0710 

LaDA0720 

LUDA0730 

LUDA0740 

UJDA0750 

LODA0760 

LODA0770 

LUDA0780 

LDDA0790 

LUDA0800 

LUDA0810 

LDDA0820 

LUDA0830 

LUDA0840 

LUDA0850 

-LUDA0860 

LUDA0870 

LUDA0880 

LDDA0890 

LDDA0900 

UIDA0910 

LUDA0920 

LODA0930 

LnDA0940 

IiODA0950 

LUDA0960 

LUDA0970 

LDDA0980 

LUDA0990 

LDDAIOOO 

LUDAIOIO 

LUDA1020 

LUDA1030 

LUDA1040 

LUDA1050 

LUDA1060 

LDDA1070 

LUDA1080 

LUDA1090 

LUDAllOO 

LUDAlllO 

LUDA1120 

LUDA1130 

LUDA1140 

LUDA1150 

LUDA1160 

LUDA1170 

LUDA1180 

LUDA1190 

LUDA1200 

LUDA1210 

LUDA1220 

LUDA1230 

LUDA1240 
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15 

20 


C 


25 


30 

35 

40 


C 


45 

50 


C 


55 


60 

65 


C 


70 


C 


75 


WI  ■  MPft 

IF  (IMX  .LT.  1)  GO  TO  20 
DO  15  K-l.lMl 

T  «  10(I,K)*L(J(K,J) 

SOM  ■  SDM-T 
NI  •  VI-t-DJVBS  (T) 

CXDBmilDB 
LOd.J)  -  SOM 
WI  »  WI-t-OABS  (SOM) 

IF  (AI  .50.  ZBRO)  AI  -  BlOft 
TEST  «  WI/AI 

IF  (TEST  .GT.  WRBEL)  WRBL  >  TEST 
GO  TO  35 

WITHCXIT  ACCORACY 
IF  (IMl  .LT.  1)  CK>  TO  35 
DO  30  K«1,IM1 

SOM  -  SOM<LO(I,K)*LO(K,J) 

CX>NTINOB 
LUd.J)  -  SOM 
cxarriNOB 
P  s  ZBRO 

COMPOTE  0(J,J)  AMD  Ld,  J)  ,  I-J+l,  . 

DO  70  I«J,R 

SOM  «  LOd.J) 

IF  (IDGT  .BQ.  0)  (30  TO  55 

WITH  ACCORACY  TEST 

AI  «  DABS (SOM) 

WI  «  ZBRO 

IF  (JMl  .LT.  1)  GO  TO  50 
DO  45  K«1,JM1 

T  -  LO(I,K)*LO(K,J) 

SOM  •  SlBl'T 
WI  •  WI+DABS(T) 

CONTINOB 
LO(I,J)  «  SOM 
WI  -  WI+DABS(SOM) 

IF  (AI  .BQ.  ZBRO)  AI  -  BI(3A 
TEST  >  WI/AI 

IF  (TEST  .GT.  WRBL)  WRBL  «  TEST 
GO  TO  65 

WITHOOT  ACCORACY  TEST 
IF  (Jm  .LT.  1)  GO  TO  65 
DO  60  K>1,JM1 

SOM  >  SOM-LO(I,K)*La(K,J) 

CONTINOB 

LO(I,J)  «  SOM 

Q  »  EQUIL(I)*DABS(SOM) 

IF  (P  .GB.  Q)  (30  TO  70 
P  >  Q 
IMAX  «  I 
CONTINOB 

TEST  FOR  ALGORITHMIC  SINGULARITY 
IF  (RN-l-P  .BQ.  RN)  GO  TO  110 
IF  (J  .BQ.  IMAX)  GO  TO  80 

INTERCHANGE  ROWS  J  AND  IMAX 

D1  «  -D1 
DO  75  K»1,N 

P  -  1U(IMAX,X:) 

LOdMAX.K)  .  LO(J,R) 

LO(J,K)  >  P 
(XRITIMOB 

‘BQOIL(IMAX)  -  BQOIL(J) 


LaDA1250 

Laua.260 

LODA1270 

L0DA1280 

L0DA1290 

L0DA1300 

LDDA1310 

LDDA1320 

L0DA1330 

LODA1340 

LODA1350 

LODA1360 

LDDA1370 

LDDA1380 

L0DA1390 

LDDA1400 

LDDA1410 

LaDA1420 

LDDA1430 

L00A1440 

LDDA1450 

,L0DA1460 

LDDA1470 

LODA1480 

LDDA1490 

LDDA1500 

LDDA1510 

L0DA1520 

L0DA1530 

LDDA1540 

LDDA1550 

L0DA1560 

LDDA1570 

LDDA1580 

LODA1590 

LDDA1600 

LDDA1610 

LDDA1620 

LDDA1630 

LDDA1640 

L0DA1650 

LODA1660 

LDDA1670 

LDDA1680 

L0DA1690 

LODA1700 

L0DA1710 

LaDA1720 

LDDA1730 

L0DA1740 

LDDA1750 

LDDA1760 

LDDA1770 

LDDA1780 

LODA1790 

LODA1800 

LDDA1810 

LDDA1820 

L0DA1830 

LODA1840 

LODA1850 

LODA1860 
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80  iFVT(j)  -  naoc 
01  •  01*10 (J,J) 

85  IF  (OABS(Ol)  .LB.  Om)  GO  TO  90 
01  •  01*SIXIR 
02  •  02+F0Cm 
GO  TO  85 

90  IP  (DABS (01)  .GB.  SIXTH)  GO  TO  95 
01  «  D1*SIZTN 
02  •  D2-F0DR 
GO  TO  90 
95  CXXnTMGB 
JPl  >  J+1 

IF  (JPl  .GT.  H)  GO  TO  105 

C  DIVIOB  BY  PIVOT  BLBMBNT  n(J,J) 

P  -  LD(J.J) 

DO  100  I>JP1,N 

LG(I,J)  -  La(I.J)/P 
100  CONTINDB 
105  CONTINDB 

C  PBRFORM  ACCURACY  TBST 

IF  (lOGT  .BQ.  0)  GO  TO  9005 
P  >  3*N-l-3 
HA  *  P*HREL 

IF  (1IA-i-10.DO**(-IOOT)  .NB.  HA)  GO  TO  9005 
IBR  «  34 
GO  TO  9000 

C  ALGC»ITKMIC  SIHQDLARITY 

110  IBR  «  129 
01  >  ZBRO 
02  •  ZBRO 
9000  CONTINDB 

C  PRINT  ERROR 

CALL  UBRTST(IBR,6HLUDATF) 

9005  RETURN 

BHD 


LDDA1870 

LTOAISSO 

LDDA1890 

LDDA1900 

LUDAIOIO 

LDDA1920 

LDDA1930 

iaDA1940 

LUDA1950 

LDDA1960 

L(IDA1970 

LDDA1980 

LDDA1990 

LDDA2000 

UJDA2010 

LDDA2020 

LDDA2030 

LaDA2040 

LaDA2050 

LDDA2060 

LODA2070 

LaOA2080 

LUOA2090 

LUDAZIOO 

LDDA2110 

LDDA2120 

L0OA2130 

LDDA2140 

LDDA2150 

LDDA2160 

LDDA2170 

LUDA2180 

LDDA2190 

LDDA2300 

LDDA2210 
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on  o  ooooooononoooooooooooooooooonooooooooooononoooo 


LlianST  RBVXSION 
PORPOSI 


USMOB 

ARGRBoanrs  a 


•  IBH/Donu 

-  JMnounr  i,  i978 

-  ILIlIXlllffIGII  PART  OP  SOLOTICai  OF  AX>B 

(POLL  STORAS8  MODI) 

-  CALL  LDSIJIP  (A,B.XPVT,N,XA,X) 

•  A  .  LQ  <THB  RSSDLT  CCMPimED  IN  THE  IMSL 

BOOTIHS  LODATP)  WmB  L  IS  A  LONBR 
TKIABnOLAR  mOSIX  WITH  0MB8  ON  IHB  MAIN 
DZAOOHAL.  a  IS  OPPSR  TRIAMSOLAR.  L  AND  U 
ARK  STORID  AS  A  SINOLB  MATRIX  A  AMD  THE 
UNIT  DIAGONAL  OP  L  IS  MOT  STCMUO).  (IHPOT) 

-  B  IS  A  VECTOR  OV  LBIRITM  M  ON  THE  RKNIT  HAND 

SIDE  OP  THE  EQONTXON  AX-B.  (IMPDT) 

•  THE  PBRMDTAnON  MAXRZX  RBTORMED  PROM  THE 

IMSL  ROOTIMB  LDDATP.  STORED  AS  AM  M  LOnTH 
VECTOR.  <IMPDT) 

-  ORDER  OP  A  AMD  MOHBSR  OP  RONS  IM  B.  (IMPOT) 

•  RON  DIMBBSION  OP  A  EXACTLY  AS  SPECIFIED  IM 

THE  OIMEMSION  STATEMEMT  IM  THE  CALLING 
PROGRAM.  (IMPOT) 

.  toe  result  X.  (OUTPUT) 

-  SIMQLB  AMD  DOUBLB/H32 

•  SIMOLB/1)36.H48,H€0 


mXCISION/MARDHARB 


REQD.  IMSL  ROUTIMES  •  MORE  REQUIRED 

NOTATION  •  IMPCMannON  OR  SPECIAL  MOXATION  AMD 

CONVEMTIONS  IS  AVAILABLE  IN  THE  MANUAL 
INTROOUCriOM  OR  THROOCRl  IMSL  ROOTTME  UHBLP 

COPYRK^  -  1978  BY  IMSL,  INC.  ALL  RIGHTS  RESERVED. 

NARRAMTY  •  IMSL  NARRAMTS  ONLY  THAT  IMSL  TESTIMO  HAS  BEBI 

APPLIED  TO  THIS  CODE.  MO  OTHER  NARRAMTY, 
EXPRESSED  OR  IMPLIED,  IS  APPLICABLE. 


SUBROUTIMB  LUBLMP  (A,B,  IPVT,M,IA,X) 

DIMERSION  A(IA,1),B(1),IPVT(1),X(1) 

DOUBLE  PRECISICW  A,B,X,SUM 

FIRST  EXECUTABLE  STATEMENT 
SOLVE  LY  »  B  FOR  Y 

DO  S  I-1,M 
5  X(I)  .  B(I) 

IN  •  0 
DO  20  I-1,M 
IP  >  IPVT(I) 

SUM  «  X(IP) 

X(IP)  -  X(I) 

IP  (IN  .BQ.  0)  GO  TO  15 
IMl  «  I-l 


LUEPOOIO 

LDV0020 

■LDBP0030 

LDBP0040 

LUEPOOSO 

LUEPOOCO 

U)EP0070 

LUEPOOSO 

LUEPOOSO 

LtnCPOlOO 

LUEPOllO 

L0EP0120 

LOEP0130 

LUEP0140 

LISPOISO 

LUEPOISO 

LUEP0170 

LUEP0180 

LUEP0190 

LUEP0200 

LDEP0310 

LDEP0220 

L0EP0230 

L0EP0240 

LDEP0250 

LDEP0260 

umro27o 

LOEP0280 

LDEP0290 

LDEP0300 

LUBP0310 

L(mP0320 

XiDEP0330 

LaEP0340 

LUEP0350 

LUEP0360 

LUEP0370 

LUBP0380 

LUEP0390 

LUEP0400 

LUEP0410 

LUBP0420 

LUBF0430 

LUEF0440 

LUEP0450 

-L0BP0460 

LUBP0470 

LUEP0480 

LUEP0490 

LUBP0500 

LUEP0510 

LUBP0S20 

LUBF0530 

L0BP0540 

LUBF0550 

LUBP0560 

LUBFOS70 

LUBP0580 

L0BP0590 

LUBP0600 

LUBP0610 

LOBF0620 


80 


c 


po  10  j-zw.na 

SDK  >  Stai-X(I.J)*X(0) 

10  cxagmoi 

GO  TO  20 

15  IF  (SOM  .NB.  O.DO)  IW  •  Z 
20  X(I)  -  SDM 

SOLVE  OX 

DO  30  ZB«1,N 
I  •  R-t-l-ZB 
IPl  .  I+l 
SOM  -  X(I) 

IP  (IPl  .or.  N)  GO  TO  30 
DO  25  J-IPl.N 

SOM  •  SDM-X(I,J)*X(J) 

25  CXWTIMUB 
30  X(l)  •  SDM/A(I,I) 

RITOim 

BHD 


LOBF0630 

UIBP0640 

LDBP0650 

LDBFOBSO 

LOBF0670 

U1BF0680 

Y  FOR  X  LDBF0S90 

LDBFOTOO 

LDBF0710 

LDBF0720 

UIBF0730 

LDBF0740 

LOBF0750 

LDBF0760 

LDBF0770 

iaBF0780 

LOBF0790 

LDBF0800 
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XMSL  ROOTUn  NMB  •  LIQTIB 


CGMPDTIR  •  IBM/DOOBLB 


lATBST  RBVISION  •  JARORRY  1,  1978 

PORPOSB  •  LINBAR  BQOATION  SOLUTION  •  BAND  STORAOB 

NODB  -  SPACB  BCCNIGMIZBR  SOLDTICW 


USACT 


•  CALL  LBQTIB  <A,N,RLC,NDC,  IA,B,M,  ZB,  IJ(»,XL, 
IBR) 


ARGDMBNTS  A 
N 

NLC 

NDC 

lA 


B 


M 

IB 


IJOB 


XL 


IBR 


-  INPUT/OOTPDT  MATRIX  OP  DIMBNSION  N  BY 

(NDC-t-NLC-t-1)  .  SBB  PARAMBTBR  IJOB. 

•  ORDBR  OP  MATRIX  A  AND  THB  NOMBBR  OP  RONS  IN 

B.  (INPUT) 

•  NDMBBR  OF  LOHBR  C0DIAG(»1ALS  IN  MATRIX  A. 

(INPUT) 

-  NOMBBR  OF  UPPBR  C0DIAG(N1ALS  IN  MATRIX  A. 

(INPUT) 

-  RON  DIMBNSICRI  OP  MATRIX  A  BXACTLY  AS 

SPBCIPIBD  IN  TUB  DIMBNSION  STATBMBNT  IN  TUB 
CALLING  PROGRAM.  (INPUT) 

-  INPUT/OOTPDT  MATRIX  OP  DIMBHSI(»I  R  BY  M. 

ON  INPUT,  B  CCXSTAINS  THB  M  RICarr-HAND  SIDBS 
OP  THB  BQDATION  AX  -  B.  OR  OUTPUT.  THB 
SOLOTICNl  MATRIX  X  RBPLACBS  B.  IP  IJOB  -  1, 

B  IS  NOT  DSHD. 

•  NOMBBR  OP  RIUHT  HAND  SIDBS  (COLOMNS  IN  B)  . 

(INPUT) 

•  RON  DIMBRSION  OP  MATRIX  B  EXACTLY  AS 

SPBCIPIBD  IN  THB  DIMBRSION  STATBMBNT  IN  THE 
CALLING  PROGRAM.  (INPDT) 

-  INPUT  OPTICm  PARAMBTBR.  IJOB  •  I  IMPLIBS  NHBR 

I  «  0,  FACTOR  THB  MATRIX  A  AND  SOLVE  THB 
BQOATION  AX  >  B.  ON  INPUT,  A  CONTAINS  THB 
COBPFICIBRT  MATRIX  OP  THB  BQOATION  AX  >  B, 
NHBRB  A  IS  ASSUMBD  TO  BB  AN  R  BY  N  BAND 
MATRIX.  A  IS  STORED  IN  BAIQ)  STORAGE  MCX3B 
AND  THEREFORE  HAS  DIMBNSION  N  BY 
(NLC-t-NUC-i-1)  .  OR  OUTPUT,  A  IS  REPLACED 
BY  THB  U  MATRIX  OP  THB  L-U  DECOMPOSITION 
OF  A  RONNISB  PBRMUTATKXT  OP  MATRIX  A.  U 
IS  STORED  IN  BAND  STORA6B  MODE. 

I  «  1,  FACTOR  THB  MATRIX  A.  A  CONTAINS  THE 
SAME  INPUT/OOTPOT  INFORMATION  AS  IF 
IJOB  -  0. 

1-2,  SOLVE  THE  EQUATION  AX  -  B.  THIS 
OPTION  IMPLIBS  THAT  LBQTIB  HAS  ALREADY 
BEEN  CALLED  USING  IJOB  -  0  OR  1  SO  THAT 
THB  MATRIX  A  HAS  ALREADY  BEEN  FACTORED. 

IN  THIS  CASE,  OUTPUT  MATRICES  A  AND  XL 
MUST  HAVB  BBBN  SAVBD  FOR  REUSE  IN  THB 
CALL  TO  LBQTIB. 

-  WORK  AREA  OF  DIMBRSION  N*(NLC-i-l)  .  THB  FIRST 

NLC^N  LOCATIONS  OF  XL  CONTAIN  COMPONENTS  OF 
THB  L  MATRIX  OP  THB  L-U  DBCCMPOSITTON  OF  A 
ROWWISE  PERMUTATION  OP  A.  THE  LAST  R 
LOCATICEIS  CONTAIN  THB  PIVOT  INDICES. 

-  ERROR  PARAMBTBR.  (OUTPUT) 


LBIBOOIO 

LB1B0020 

■LE1B0030 

LB1B0040 

LB1B0050 

LB1B0060 

LB1B0070 

LB1B0080 

LB1B0090 

LBIBOIOO 

UtlBOllO 

LB1B0120 

LB1B0130 

LB1B0140 

LBIBOISO 

LB1B0160 

LB1B0170 

LB1B0180 

LB1B0190 

LB1B0200 

LB1B0210 

LB1B0220 

LB1B0230 

LB1B0240 

LB1B0250 

LB1B0260 

LB1B0270 

LB1B0280 

LB1B0290 

LB1B0300 

LB1B0310 

LB1B0320 

LB1B0330 

LB1B0340 

LB1B0350 

LB1B0380 

LB1B0370 

LB1B0380 

LB1B0390 

LB1B0400 

LB1B0410 

LB1B0420 

LB1B0430 

LB1B0440 

LB1B0450 

LB1B0440 

LB1B0470 

LE1B0480 

LB1B0490 

LB1B0500 

LB1B0510 

LB1B0520 

LB1B0S30 

LB1B0540 

LEIBOSSO 

LB1B0560 

LE1B0570 

LE1B0580 

LB1B0590 

LB1B0600 

LBlBOeiO 

LB1B0620 
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nnn  n  n  nnnnnnnnnnnnnnnnnnnnnn 


PRBCISION/HARDHARB 


TSRMIKAL  BRftOR 

XIR  -  139  UDICATBS  THAT  MATRIX  A  IS 
AI^RITHMICALLY  SIMQDLAR.  (SKI  THE 
CHAPTIR  1.  PRXLDDK)  . 

SINGLB  AMD  DCX1BLB/H32 
SINQLR/H3€,H48,H€0 


RBQD.  IMSL  ROOTTRBS  >  DBRTST.OQKTIO 


NOTATION 


CX>PYRI6HT 


HARRANTT 


INFORMATIW  ON  SPECIAL  NOTATION  AND 
CONVENTIONS  IS  AVAILABLE  IN  THE  MANUAL 
INTRQDOCTION  OR  THROUGH  IMSL  ROUTINE  UHELP 

1978  BY  IMSL,  INC.  ALL  RIGHTS  RESERVED. 

IMSL  WARRANTS  ONLY  THAT  IMSL  TESTING  HAS  BEEN 
APPLIED  TO  THIS  C(X)E.  HO  OTHER  WARRANTY, 
EXPRESSED  OR  IMPLIED,  IS  APPLICABLE. 


SUBROUTINE  LECJTIB  (A,N,NLC,NUC,  IA,B,M,  IB,  IJOB.XL,  lER) 


DIMENSION 
DOUBLE  PRECISION 
DATA 


A(IA,1)  ,XL(N,1)  ,B(IB,1) 
A,XL,B.P,Q,ZERO,ONB,RN 
ZERO/O .  DO/ ,  CNIE/1 .  ODO/ 

FIRST  EXECUTABLE  STATEMENT 


lER  «  0 
JBEG  »  NLC-fl 
NLCl  •  JBEG 

IF  (IJCS  .EQ.  2)  GO  TO  80 
RN  «  N 


I  «  1 

NC  «  JBEG-t-NUC 
NN  ■  NC 
JEND  «  NC 
IF  (N  .EQ.  1  . 
5  K  «  1 
P  -  ZERO 
DO  10  J  >  JBK( 


RESTRUCTURE  THE  MATRIX 

FIND  RECIPROCAL  OF  THE  LARGEST 

ABSOLUTE  VALUE  IN  ROW  I 


OR.  NLC  .EQ.  0)  GO  TO  25 


DO  10  J  -  JBEG, JEND 
A(I,K)  ■  A(I,J) 

Q  -  DABS(A(I,K)) 

IF  (Q  .GT.  P)  P  .  Q 
K  «  K-t-1 
10  CONTINDE 

IF  (P  .EQ.  ZERO)  GO  TO  135 
XL(I,NLC1)  -  ONE/P 
IF  (K  .GT.  NC)  GO  TO  20 
DO  15  J  «  K,RC 
A(I,J)  -  ZERO 
15  CONTINUE 
20  I  >  I<fl 

JBEG  -  JBKG-1 

IF  (JKND-JBKG  .EQ.  N)  JEND  »  JKRD-1 
IF  (I  .LE.  NLC)  GO  TO  5 
JBEG  «  I 
NH  *  JEND 
25  JKto  >  N-NUC 


LK1B0630 

LK1B0640 

LX1B0650 

LK1B0660 

LE1B0670 

LK1B0680 

LK1B0690 

LE1B0700 

LE1B0710 

LK1B0720 

LK1B0730 

LE1B0740 

LE1B0750 

LE1B0760 

LK1B0770 

LE1B0780 

LE1B0790 

LE1B0800 

LE1B0810 

LE1B0820 

-LB1B0830 

LE1B0840 

LE1B0850 

LE1B0860 

LK1B0870 

LE1B0880 

LK1B0890 

LE1B0900 

LK1B0910 

LK1B0920 

LK1B0930 

LK1B0940 

LK1B0950 

LE1B0960 

LE1B0970 

LK1B0980 

LK1B0990 

LKIBIOOO 

LKIBIOIO 

LK1B1020 

LE1B1030 

LE1B1040 

LE1B1050 

LK1B1060 

LE1B1070 

LE1B1080 

LE1B1090 

LEIBIIOO 

LElBlllO 

LE1B1120 

LE1B1130 

I.E1B1140 

LE1B1150 

LE1B1160 

LE1B1170 

LE1B1180 

LE1B1190 

LE1B1200 

LE1B1210 

LE1B1220 

LE1B1230 

LE1B1240 
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DO  40  I  >  JBBQ.N 

LB1B1250 

P  -  ZBRO 

IJI1B1260 

DO  30  J  >  1,101 

LB1B1270 

Q  »  DABS(A(I,J)) 

LB1B1280 

IP  (Q  .GT.  P)  P  .  Q 

Lb1B1290 

30 

cxorrniDB 

LB1B1300 

IF  (P  .BQ.  ZBRO)  GO  TO  135 

LB1B1310 

XLd.NLCl)  •  ONB/P 

LB1B1320 

IF  (I  .BQ.  JBRD)  GO  TO  37 

LB1B1330 

IF  (I  .LT.  JBND)  GO  TO  40 

LB1B1340 

K  «  ini-»-i 

LB1B1350 

DO  35  J  -  K,NC 

LB1B1360 

A(I,J)  »  ZBRO 

LB1B1370 

35 

CONTIMUB 

LB1B1380 

37 

NN  »  NN-l 

LB1B1390 

40 

CORTTHDB 

LB1B1400 

L  >  NLC 

LB1B1410 

L-U  DECOMPOSITION 

LB1B1420 

DO  75  K  .  1,H 

LB1B1430 

P  «  DABS(A(K.l) )*XL(K,NLC1) 

LB1B1440 

I  -  K 

LB1B1450 

IP  (L  .LT.  N)  L  «  Irfl 

LB1B1460 

K1  >  K+l 

LB1B1470 

IF  (K1  .GT.  L)  GO  TO  50 

LB1B1480 

DO  45  J  «  K1,L 

LB1B1490 

Q  >  DABS(A(J,1) )*XL(J,NLC1) 

LB1B1500 

IF  (Q  .LB.  P)  GO  TO  45 

LB1B1510 

P  >  Q 

LB1B1520 

I  >  J 

LB1B1530 

45 

CORTINDB 

LB1B1540 

50 

XL(I,NLC1)  •  XL(K,NLC1) 

LB1B15S0 

XL(K,NLC1)  •  I 

LB1B1560 

SINGOLARXTT 

FOOND 

LB1B1570 

Q  -  RN-fP 

LB1B1580 

IP  (Q  .BQ.  RH)  GO  TO  135 

LB1B1590 

INTBRCHAIIOB 

ROWS  I  AND  K 

LB1B1600 

IF  (K  .BQ.  I)  GO  TO  60 

LB1B1610 

DO  55  J  -  1,NC 

LB1B1620 

P  «  A(K,J) 

LB1B1630 

A(K,J)  .  A(I,J) 

LB1B1640 

A(I,J)  .  P 

LB1B1650 

55 

CONTINUB 

LB1B1660 

60 

IF  (K1  .GT.  L)  GO  TO  75 

LB1B1670 

DO  70  I  «  K1,L 

LB1B1680 

P  -  A(I,1)/A{K,1) 

LB1B1690 

IK  -  I-K 

LB1B1700 

XL(K1,IK)  «  P 

LE1B1710 

DO  65  J  >  2,NC 

LB1B1720 

A{I.J-1)  -  A(I,J) -P*A(K,J) 

LB1B1730 

65 

CCBfriNDE 

LB1B1740 

A(I,RC)  *  ZBRO 

LB1B1750 

70 

CONTINDB 

LB1B1760 

75 

CONTINDB 

LB1B1770 

IF  (IJOB  .BQ.  1)  GO  TO  9005 

LE1B1780 

FORNARD  SOBSTITDTION 

LB1B1790 

80 

L  «  NLC 

LB1B1800 

DO  105  K  «  l.N 

LE1B1810 

I  >  XL(K,NLC1) 

LB1B1820 

IF  (I  .BQ.  K)  GO  TO  90 

LB1B1830 

DO  85  J  -  1,M 

LB1B1840 

P  >  B(K,J) 

LB1B1850 

B(K,J)  «  B(I,J) 

LB1B1860 
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c 


B(I.J)  «  P 

LB1B1870 

85 

COMTINDB 

LB1B1880 

90 

ZP  (1.  .LT.  N)  L  >  L+l 

LB1B1890 

K1  •  K-fl 

LB1B1900 

IP  (Kl  .GT.  L)  GO  TO  105 

LE1B1910 

DO  100  I  -  K1,L 

LE1B1920 

IK  ■  I-K 

LE1B1930 

P  >  XL(K1,IK) 

LB1B1940 

DO  95  J  «  l.M 

LB1B1950 

B(I,J)  -  B(I,J) -P*B(K.J) 

LB1B1960 

95 

CCBmNDB 

LB1B1970 

100 

CONTIRDB 

LB1B1980 

105 

CONTINUE 

LB1B1990 

BACKMARD  SUBSTITUTION 

LB1B2000 

JBBG  «  NUC+NLC 

LB1B2010 

DO  125  J  ■  l.M 

LB1B2020 

L  ■  1 

LB1B2030 

Kl  «  N-t-1 

LB1B2040 

DO  120  I  >  1,N 

LE1B2050 

K  >  Kl-I 

LE1B2060 

P  >  B(K,J) 

LE1B2070 

IP  (L  .BQ.  1)  GO  TO  115 

LB1B2080 

DO  110  KK  -  2,L 

LB1B2090 

IK  >  KK-t-K 

LB1B2100 

P  -  P-A(K,KK)*B (IK-1, J) 

LB1B2110 

110 

CONTINUE 

LE1B2120 

115 

B(K,J)  -  P/A(K,1) 

LB1B2130 

IP  (L  .LB.  JBEQ)  L  •  lri-1 

LB1B2140 

120 

CONTINUE 

LB1B2150 

125 

CONTINUE 

LB1B21C0 

GO  TO  9005 

LB1B2170 

135 

lER  >  129 

LB1B2180 

9000 

CONTINUE 

LB1B2190 

CALL  UBRTST(IBR,6HLEQT1B) 

LB1B2200 

9005 

RETURN 

LB1B2210 

END 

LB1B2220 
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msL  RODTnn  miB  •  obstst 


COHraTBR 
latest  Rsvisicai 


PURPOSE 


USAGE 


ARGUMENTS 


NAME 


•  ZBM/SINGLE 

-  JUNE  1,  I9S2 

-  PRINT  A  MESSAGE  REFLECTING  AH  ERROR  CONDITION 

-  CALL  UERTST  (ISR.NAME) 

•  ERROR  PARAMETER.  (INPUT) 

lER  -  I+J  NHSRE 

I  -  128  IMPLIES  TERMINAL  ERROR  MESSAGE, 

I  «  64  IMPLIES  WARNING  WITH  FIX  MESSAGE, 

I  >  32  IMPLIES  WARNING  MESSAGE. 

J  -  ERROR  CODE  RELEVANT  TO  CALLING 
ROUTINE. 

-  A  CHARACTER  STRING  OP  LENGTH  SIX  PROVIDING 

THE  NAME  OF  THE  CALLING  ROUTINE.  (INPUT) 


PRECISION/KARDWARE  -  SINGLE/ALL 
RBQD.  IMSL  ROUTINES  -  UGETIO,USPKD 

NOTATION  -  INFORMATION  ON  SPECIAL  NOTATION  AND 

CONVENTIONS  IS  AVAILABLE  IN  THE  MANUAL 
INTRODUCTION  OR  THROUGH  IMSL  ROUTINE  UHELP 

REMARKS  THE  ERROR  MESSAGE  PRODUCED  BY  UERTST  IS  WRITTEN 

TO  THE  STANDARD  OUTPUT  UNIT.  THE  OUTPUT  UNIT 
NUMBER  CAN  BE  DETERMINED  BY  CALLING  (A3BTIO  AS 
FOLLOWS. .  CALL  UGBTIO(l,NIN,NODT) . 

THE  OUTPUT  UNIT  NUMBER  CAN  BE  CHANGED  BY  CALLING 
UGETIO  AS  FOLLOWS . . 

NIN  «  0 

NOUT  •  NEW  OUTPUT  UNIT  NUMBER 
CALL  UGETIO (3, NIN, NOUT) 

SEE  THE  UGETIO  DOCUMENT  FOR  MORE  DETAILS. 

COPYRIGHT  -  1982  BY  IMSL,  INC.  ALL  RIGHTS  RESERVED. 

WARRANTY  -  IMSL  WARRANTS  ONLY  THAT  IMSL  TESTING  HAS  BEEN 

APPLIED  TO  THIS  COOE.  NO  OTHER  WARRANTY, 
EXPRESSED  OR  IMPLIED,  IS  APPLICABLE. 


SUBROUTINE  UERTST  (IER,NAME) 

SPECIFICATIONS  FOR  ARGUMENTS 

INTEGER  lER 

INTEGER  NAMB(l) 

SPECIFICATIONS  FOR  LOCAL  VARIABLES 
INTEGER  I,IEQ,IEQDF,IOUNIT,LBVBL,LEVOLD,NAMEQ(6) , 

>  NAMSET(6)  ,NAMUPK(6)  ,NIN,NMTB 

DATA  NAMSET/1HU,1HE,1KR,1HS,1HE,1HT/ 

DATA  NAMEQ/6*1H  / 

DATA  LBVEL/4/ , IBQDF/0/, IEQ/1H«/ 

UNPACK  NAME  INTO  NAMUPK 
FIRST  EXECUTABLE  STATEMENT 
CUVLL  USPKD  (NAME,  6, NAMUPK, NMTB) 


UERTOOIO 
UBRT0020 
-UBRT0030 
UBRT0040 
UERTOOSO 
UBRT0060 
UBRT0070 
UERTOOSO 
UERT0090 
UERTOlOO 
UERTOllO 
UERT0120 
UERT0130 
UBRT0140 
UERT0150 
UERT0160 
UERT0170 
UERT0180 
UERT0190 
UERT0200 
UERT0210 
UERT0220 
UERT0230 
UERT0240 
UERT0250 
UERT0260 
UERT0270 
UERT0280 
UERT02S3 
UERTOSP*^ 
UERT0.31  . 
UERT0320 
UERr0330 
UERT0340 
UERT0350 
UERT0360 
UERT0370 
UERT0380 
UERT0390 
UERT0400 
UBRT0410 
UBRT0420 
UERT0430 
UBRT0440 
UBRT0450 
UERT0460 
UERT0470 
-UBRT0480 
UBRT0490 
UERTOSOO 
UBRT0510 
UBRT0520 
UBRT0530 
UERT0540 
UBRT05S0 
UBRT0560 
DERT0570 
UBRT0580 
UBRT0590 
UBRT0600 
UBRT0610 
UERT0620 
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onno 


C 


GET  OOTPtTT  UNIT  NOMBBR 


CALL  GQBTIO(1.NIN,IOONIT) 

C  CHECK  lER 

IF  (ZRR.GT.999)  GO  TO  25 

IF  (IER.LT.>32)  GO  TO  55 

IF  (IER.LB.128)  GO  TO  5 

IF  (LBVEL.lt. 1)  GO  TO  30 

C  PRINT  TERldNAL  MESSAGE 

IF  (IBQDF.BQ.l)  WRITE (lOONIT. 35)  IBR, NAMBQ, IBQ, NAMDPK 
IF  (IBQDF.BQ.O)  WRITE (lOONIT, 35)  IBR.MAMDPK 
GO  TO  30 

5  IF  (IER.LB.64)  GO  TO  10 
IF  (LBVBL.lt. 2)  GO  TO  30 

C  PRINT  WARNING  WITH  FIX  MESSAGE 

IF  (IBQDF.BQ.l)  WRITE (lOONIT, 40)  IBR.NRMEQ, IBQ.NAMDPK 
IF  (IBQDF.BQ.O)  WRITE (lODNIT. 40)  IBR.NAMDPK 
GO  TO  30 

10  IF  (IBR.LB.32)  GO  TO  15 

C  PRINT  WARNING  MESSAGE 

IF  (LEVEL.lt. 3)  GO  TO  30 

IF  (IBQDF.BQ.l)  WRITE (lOONIT, 45)  lER.NAMBQ, IEQ,HAMDPK 
IF  (IBQDF.BQ.O)  WRITE (lOONIT, 45)  IBR.NAMOPK 
GO  TO  30 
15  CONTINUE 

C  CHECK  FOR  DBRSBT  CJOL 

DO  20  I>1,6 

IF  (NAMUPK(I) .NE.HAMSBT(I) )  GO  TO  25 
20  CONTINUE 

LBVOLD  >  LEVEL 
LEVEL  >  IBR 
lER  a  LEVOLD 

IF  (LEVEL.lt. 0)  LEVEL  •  4 
IF  (LEVEL. GT. 4)  LEVEL  «  4 
GO  TO  30 
25  CONTINUE 

IF  (LEVEL.lt. 4)  GO  TO  30 

C  PRINT  NON-DBFINED  MESSAGE 

IF  (IBQDF.BQ.l)  WRITE (lODNIT, 50)  IBR.NAMBQ, lEQ.NAMUPK 
IF  (IB^F.EQ.O)  WRITE (lOUNIT, 50)  IER,NAMUPK 
30  lEQDF  «  0 
RETURN 

35  FORMAT  (19H  ***  TERMINAL  ERROR,  1  OX,  7H  (lER  =  ,  13  , 

1  20H)  FRCXf  IMSL  ROUTINE  ,€A1,A1,6A1) 

40  FORMAT  (27H  ***  WARNING  WITH  FIX  ERROR,  2X,  7H  (lER  »  ,13, 

1  20H)  FROM  IMSL  ROUTINE  ,6A1,A1,6A1) 

45  FORMATdSH  ***  WARNING  ERROR,  1 IX, 7H(IER  »  ,13, 

1  20H)  FRC»1  IMSL  ROUTINE  ,6A1,A1,6A1) 

50  FORMAT (20H  ***  UNDEFINED  ERROR, 9X, 7H(IER  »  ,15, 

1  20H)  FROM  IMSL  ROUTINE  ,fiAl,Al,6Al) 

SAVE  P  FOR  P  -  R  CASE 
P  IS  THE  PAGE  NAHUPK 
R  IS  THE  ROUTINE  NAMUPK 

55  lEQDF  s  1 
DO  60  I«l,6 

60  NAMEQ(I)  >  NAMUPK(I) 

65  RETURN 
END 


UBBT0630 

UBRT0640 

UBirr0650 

'TERT0660 

UBRT0670 

UBRT0680 

UBRT0690 

DERT0700 

DBRT0710 

UBRT0720 

UBRT0730 

UBRT0740 

UERT0750 

UERT0760 

DBRT0770 

DBRT0780 

UBRT0790 

UBRT0800 

UERT0810 

UERT0820 

UERT0830 

UERT0840 

UERT0850 

UERT0860 

UBRT0870 

UERT0880 

UBRT0890 

UBRT0900 

UERT0910 

UBRT0920 

UBRT0930 

UBRT0940 

UBRT0950 

UERT0960 

UBRT0970 

UBRT0980 

UERT0990 

UBRTIOOO 

UERTlOlO 

UERT1020 

UERT1030 

UERT1040 

UERT1050 

UERT1060 

UERT1070 

UERT1080 

UERT1090 

UERTllOO 

UERTlllO 

UERT1120 

UERT1130 

UERT1140 

UERT1150 

UBRT1160 

DERT1170 

UERT1180 

UERT1190 

UERT1200 
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c 

c 

c- 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 
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IMSL  RODTIMB  HAMS  -  OQBTIO 


CXMraTBR 
LATEST  REVISION 
PURPOSE 


•  IBM/SINGLE 

•  JUNE  1,  1981 

-  TO  RETRIEVE  CURRENT  VALUES  AND  TO  SET  NEW 
VALUES  FOR  INPUT  AND  OUTPUT  UNIT 
IDENTIFIERS. 


USAGE 


CALL  UGETIOdOPT.NIN.NOUT) 


arguments  IOPT 


NIN 

NOUT 


OPTION  PARAMETER.  (INPUT) 

IF  IOPT-1,  THE  CURRENT  INPUT  AND  OUTPUT 
UNIT  IDENTIFIER  VALDES  ARE  RETURNED  IN  NIN 
AMD  NOUT,  RESPECTIVELY. 

IF  IOPT«2,  THE  INTERNAL  VALUE  OF  NIN  IS 
RESET  FOR  SUBSEQUENT  USE. 

IF  IOPT«3,  THE  INTERNAL  VALUE  OF  NOUT  IS 
RESET  FOR  SUBSEQUENT  USE. 

INPUT  UNIT  IDENTIFIER. 

OUTPUT  IF  IOPT-1,  INPUT  IF  IOPT- 2 . 

OUTPUT  UNIT  IDENTIFIER. 

OUTPUT  IF  IOPT-1,  INPUT  IF  IOPT- 3 . 


PRECISION/KARDMARE  •  SINGLE/ALL 


RE(^.  nCSL  ROUTINES  •  NONE  REQUIRED 

NOTATION  •  INFORMATION  ON  SPECIAL  NOTATION  AMD 

CONVENTIONS  IS  AVAILABLE  IN  THE  MANUAL 
INTRODUCTION  OR  THROUGH  IMSL  ROUTINE  UHELP 

REMARKS  EACH  IMSL  ROUTINE  THAT  PERFORMS  INPUT  AMD/OR  OUTPUT 

OPERATIONS  CALLS  UGETIO  TO  OBTAIN  THE  CURRENT  TT  ^ 
IDENTIFIER  VALUES .  IF  UGETIO  IS  CALLED  WITH  lOF  OR 
IOPT-3 ,  NEW  UNIT  IDENTIFIER  VALUES  ARB  BSTABLISh..^  . 
SUBSEQUENT  INFDT/OUTPUT  IS  PERFORMED  ON  THE  NEW  UNITS . 


COPYRIGHT 


1978  BY  IMSL,  INC.  ALL  RIGHTS  RESERVED. 


UGETOOIO 

UGBT0020 

-UCaT0030 

UGBT0040 

UGETOOSO 

OGETOOeO 

UGBT0070 

UGBT0080 

UGETOOSO 

UGETOlOO 

UGETOllO 

UGET0120 

UGBT0130 

UGET0140 

UGET0150 

0(^0160 

GGET0170 

UGET0180 

UGET0190 

UGET0200 

UGET0210 

ucarroiio 

DGET0230 

UGET0240 

UQET0250 

UGET0260 

UQET0270 

UGET0280 

UGET0290 

06ET0300 

OQET0310 

UGET0320 

tlSET0330 

UGET0340 

UQET03S0 

U6ET0360 

IK3ET0370 

UGET0380 

UGET0390 

UGET0400 

UGET0410 

UGET0420 

UGET0430 


WARRANTY  ■  IMSL  WARRANTS  ONLY  THAT  IMSL  TESTING  HAS  BEEN  (K3ET0440 

APPLIED  TO  THIS  CODE.  NO  OTHER  WARRANTY,  UGET0450 

EXPRESSED  OR  IMPLIED,  IS  APPLICABLE.  UGBT0460 

UGET0470 

. . UGET0480 

UGBT0490 

SUBROUTINE  UGETIO (IOPT, NIN, NOUT)  UGETOSOO 

SPECIFICATIONS  FOR  AR(R]MENTS  UGETOSIO 


INTEGER 

IOPT,  NIN,  NOUT 

UGBT0520 

SPECIFICATIONS  FOR  LOCAL  VARIABLES 

UGET0530 

INTEGER 

RIND, NOUTD 

UGET0540 

DATA 

NlND/5/,NOOTD/6/ 

DGET0550 

FIRST  EXECUTABLE  STATEMENT 

UGET0560 

IF  (I0PT.EQ.3) 

GO 

TO  10 

UGET0570 

IF  (I0PT.EQ.2) 

GO 

TO  5 

UGET0580 

IF  (lOPT.NE.l) 

GO 

TO  9005 

UGET0590 

NIN  -  NIMD 

UGETOSOO 

NOUT  -  NOUTD 

UGETOSIO 

GO  TO  9005 

UGET0S20 
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o  in 


5  NIMD  -  NIN 
GO  TO  9005 
NOUTD  ■  NODT 
RnURN 
BBOD 


osrro€3o 

aQBT0C40 

oQrro650 

IX»T0660 

OQBTOSTO 
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on  o  n  n  no  no  o  o  no  o  n  no  o  o  o  n  n  n  o  n  o  o  n  o  o  o  o  o  o  o  o  o  o  o  no  o 


msL  nooniB  hams  •  dspxd 


CXHODTSR  -  IBN/SIHOXiB 

lATBST  SBVISZOH  -  HOVBMBBR  1,  1984 

FCRPOSB  -  NDCLBDS  CALLBD  BY  ZHSL  ROaTINBS  THAT  HAVE 

CRAKACTBR  STRIHO  AnOOlIBBTS 

nSAQB  -  CALL  nSPKD  ( PACKED, HCKARS.UHPAKD.NCIMIB) 

ARGDUBNTS  PACKED  -  CHARACTER  STRIHO  TO  BE  DHPACKED.  (ZNPOT) 
HCHARS  -  LEHOIH  OF  PACKED.  (IHPDT)  SEE  REMARKS. 
OHPAKD  •  IirrBOSR  ARRAY  TO  RECEIVS  THE  DHPACKED 
REPRESERTATIOn  OF  THE  STRIHO.  (OOTPOT) 

.  HCHMIB  -  HCHARS  MIHOS  TRAILIHO  BLAHKS.  (ODTPDT) 

PRECISIQH/HARDIIARB  •  SIHOLB/ALL 


RB0D.  IMSL  ROC7TIHBS  -  HOHE 

REMARKS  1 .  DSPKD  DHPACKS  A  CHARACTER  STRIHO  IHTO  AH  IHTBOBR  ARRAY 
IH  (Al)  FORMAT. 

2.  DP  TO  129  CHARACTERS  MAY  BE  DSBD.  AHY  IH  EXCESS  OF 
THAT  ARE  IGHORED. 


COPYRIOHT  >  1984  BY  IMSL,  IHC.  ALL  RKHTTS  RESERVED. 

HARRANTY  •  IMSL  KARRAHTS  ONLY  THAT  IMSL  TBSTIHO  HAS  BBBH 

APPLIED  TO  THIS  CODE.  HO  OTHER  KARRAHTY, 
EXPRESSED  OR  IMPLIED,  IS  APPLICABLE. 


SDBROOTIHB  DSPKD  (PACKED, HCHARS, DHPAKD,I«CHMTB) 

SPBCIFZCATIOHS  FOR  ARGDMBHTS 
IHTBOBR  HC, HCHARS, HCHMTB 


LOGICAL*!  DHPAKD(l) , PACKED (1) ,LBYTE,LBLAHK 

IHTBQBR*2  IBYTE,IBLAHK 

SgUIVALEHCE  (LBYTE , IBYTB) 

DATA  LBLAHK  /IH  / 

DATA  IBYTE  /IH  / 

DATA  IBLAHK  /IH  / 

IHITZALIZE  HCHMTB 


HCHMTB  >  0 


IF (HCHARS. LE.O)  RBTDRH 


RBTDRN  IF  HCHARS  IS  LE  ZERO 


SET  HC.HDMBBR  OF  CHARS  TO  BE  DECODED 

HC  -  MIHO  (129, HCHARS) 

HNORDS  «  HC*4 
J  ■  1 

DO  110  I  -  l,HHORDS,4 
DHPAKD(I)  -  PACKED (J) 

DHPAKD(I+1)  -  LBLAHK 
DHPAKD(I-t-2)  «  LBLAHK 
DHPAKD(I+3)  -  LBLAHK 
110  J  >  J-t-1 


DO’ 200  H  >  l,HHORDS,4 


CHECK  DHPAKD  ARRAY  AHD  SET  HCHMTB 
BASED  OH  TRAILIHO  BLAHKS  FCXJHD 


DSPKOOlO 

DSPK0020 

-Dsnmoso 

DSPK0040 

DSPK0050 

DSPK0060 

DSPK0070 

DSPKQ080 

DSPK0090 

DSPKOlOO 

DSPKOllO 

DSPK0120 

DSPK0130 

OSPK0140 

DSPKOISO 

OSPK0160 

DSPK0170 

DSPKOISO 

DSPK0190 

DSPK0200 

DSPK0210 

DSPK0220 

DSPK0230 

DSPK0240 

D8PK0250 

DSPK0260 

DSPK0270 

DSPK0280 

DSPK0290 

DSPK0300 

DSPK0310 

DSPK0320 

DSPK0330 

DSPK0340 

>DSnC0350 

DSPK0360 

DSPK0370 

DSPK0380 

DSPK0390 

DSPK0400 

OSPK0410 

DSPK0420 

DSPK0430 

DSPK0440 

DSPK0450 

DSPK0460 

DSPK0470 

DSPK0480 

DSPK0490 

DSPK0500 

DSPK0510 

DSPK0520 

DSPK0530 

DSPK0540 

DSPK0550 

DSPK0560 

DSPK0570 

DSPK0580 

DSPK0590 

DSPK0600 

DSPK0610 

DSPK0620 


im  •  imOilDS  -  N  -  2  iraPK0630 

zam  -  OHPAXD(IDI)  nSPK0640 

IFdBTTB  .NB.  IBLANX)  GO  TO  210  USPK0650 

200  CXWTIRDB  nSPK0660 

MN  •  0  aSPX0670 

210  NCKMIB  -  (RK  -t-  3)  /  4  USPK0680 

RSTOPR  USPK0690 

mo  C«PK0700 
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